Accepted by ApJS 



Three Dimensional Magneto Hydrodynamical Simulations of 
Gravitational Collapse of a 15M Q Star 

Takami Kuroda and Hideyuki Umeda 

Department of Astronomy, School of Science, University of Tokyo, Bunkyo-ku, Tokyo, 113-0033, 
Japan; kuroda@astron.s.u-tokyo. ac.jp, umeda© astron.s.u-tokyo. ac.jp, 

ABSTRACT 

We introduce our newly developed two different, three dimensional magneto hydro- 
dynamical codes in detail. One of our codes is written in the Newtonian limit (NMHD) 
and the other is in the fully general relativistic code (GRMHD). Both codes employ 
adaptive mesh refinement and, in GRMHD, the metric is evolved with the "Baumgarte- 
Shapiro-Shibata-Nakamura" formalism known as the most stable method at present. 
We did several test problems and as for the first practical test, we calculated gravita- 
tional collapse of a 15M & star. Main features found from our calculations are; (1) High 
velocity bipolar outflow is driven from the proto-neutronstar and moves through along 
the rotational axis in strongly magnetized models; (2) A one-armed spiral structure 
appears which is originated from the low-|T/W| instability; (3) By comparing GRMHD 
and NMHD models, the maximum density increases about ~ 30% in GRMHD models 
due to the stronger gravitational effect. These features agree very well with previous 
studies and our codes are thus reliable to numerical simulation of gravitational collapse 
of massive stars. 

Subject headings: Methods: numerical — MHD — stars: magnetars — supernovae: 
general 



1. Introduction 

There are lots of works searching for the explosion mechanisms of core-collapse supernovae 
(CCSNe), however we have still not obtained any conclusive results in these decades. Recent 
works, both observational and theore tical ones, sh o w sev eral indications that their explosions are 



commonly aspherical. For instance, iMaeda et al.l (|2008l ) obtained late-time spectra for a lot of 



CCSNe and showed that the explosion morphologies of stars without H envelope are close to bipo- 
lar configu rations. Furthermore , non-axisymmetric explosion is found from the observation of SN 



2005bf by iTanaka et al.l (|2009bl ). In their report, they presented an optical spectropolarimetric 



observation of Type lb supernova 2005bf and claimed that SN 2005bf can be explained as unipolar 
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explosion and also the direction of launched unipolar blob is tilted from the symmetric axis. There- 
fore, the asphericities might be key ingredients to understand the explosion mechanisms, especially 
for a subset class of CCSNe such as Type lb SNe. These asphericities found from the observa- 
tions are thought to be products of hydrodynamical instabilities occurring in the vicinity of the 
proto neutronstars (PNSs). From the previous theoretical/numerical works, it is widely known that 
there are many types of hydrodynamica l instabilities which would occ ur during CCSNe, s uch as: 
the Ledoux convection (IKeil et al.lll996l ): vortical-acoustic instability (jBlondin et al.l 120031): mag- 
neto hydrodynamical instabilities, e.g., Magneto Rot ational Instability (MRI) (jBalbus &: Hawlev 
199ll : lAkivama et "all I2OO3I : lobergaulinger et al.ll2009l) or Kelvin Helmholtz instability: rot ational 



instabilities, e. g., dynamical bar- m ode instability (jRampp et al.ll l998; 



lar instability (Imamura et al.ll2003l ). low-|T/W| instability (jShibata et al 



Shibata et al 



2003 



2003), secu- 



Watts et al. 



2005 



Ott et al.ll2005l ). Which of these instabilities would occur depends on progenitor mass, the rota- 
tional/magnetic field velocity configuration, the interaction between matters and neutrinos, etc. 

Among these mechanisms, rotational instabilities are common byproducts of relatively fast- 
spinning progenitors and their subsequent collapses. If ratio of rotational to gravitational potential 
energy (3 = \T/W\, at core bounce, exceeds 0dy n ~ 0.27 then the dynamical bar-mode instability 



appears (jRampp et al 



1998 



Shibata et al.ll2003l ). Simple estimate, in which we assume the angular 



momentum is almost conserved during core collapse, gives us a rough lower limit of the central 
angular velocity f2 c , at pre-collapse stage which exceeds ~ 10 rad/s. Such rotati onal speed is at 
least ~ 10 times faster than the results of recent stellar evolutional calculations (jYoon &; Langer 



20051 ). Even if progenitor does not spin so rapidly at the beginning, the secular instability may be 
appeared if /3 > (3 sec ~ 0.14 and also if some dissipative m echanisms exist such as the viscosity, the 
neutrino radiation or the gravitational radiation reaction (jlmamura et al.ll2003l ). Recently, another 
i nteresting rotationa l instability is reported to occur in CCSNe which is the low-|T/W| instability 
(jShibata et al.1 120031 : IWatts et al.ll2005l : lOtt et al.ll2005l ). This is the resonance instability across 
the corotation point inside the differentially rotating area and occurs with more reasonable value 
/3 ~ 0.01 or with relatively slow initial spin rate S7 C ~ 1 rad/s. Such a slow rotation is more realistic 
compared to aforementioned two mechanisms. 

These rotational instabilities are intrinsically three dimensional, non-axisymmetric phenom- 
ena. In the ideal MHD limit, which is a reasonable assumption in the CCSNe context, the non 
axisymmetric motion always convert the toroidal magnetic field into the poloidal component by 
dragging the toroidal magnetic field. Then, if the differential rotation exists, the converted poloidal 
magnetic field is again converted back into the toroidal ones. Such a closed cycle never takes 
place in the axisymmetric mo tion and may play importa nt roles to amplify magnetic field via, 
e.g., the dynamo mechanism ([Thompson Duncan! Il993l ) or the MRI. Strongly a mplified mag- 
netic field (> 10 15 ~ 16 G) launches high velocity outflow along the rotational axis ( Mikami et al 



2008) due to the magneto-spring or the magneto-centrifugal effects (IWheeler et al.ll2002l ) and may 
l eave highly magnetized (~ 10 15_16 G at surface) neutronstar which is a so-called "magnetar" 



(jDuncan &; Thompson! Il992l ). Another feature of the magneto-rotational explosion is, through 
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these magneto-rotational effects, the explosion morphology becomes highly aspherical. Since the 
explosion morphology shows stronger asphericity in ca se of energetic exp losion such as Hypernova 
(HN) or SN associated with gamma-ray burst (GRB) ([Maeda et al.ll2008l ). the magneto rotational 
effects are considered to be more important to such hyper energetic CCSNe than the normal ones. 

As just described, the non axisymmetric effects may play important roles and should be exam- 
ined by three dimensional numerical simulations. However, up to the pres ent date, there are only a 
few n u merical works about CCS Ne with three dimensional MHD (see, e.g.. lScheidegger et al.l (|2008l . 
20091 ) ; iMikami et al.l (|2008l )). In IScheidegger et al.l (|2009l ). they calculated a number of numerical 
models aiming for the gravitational wave signature during core collapse, including two types of 
realistic EOSs, various magnetic-rotational configurations and a neutrino parametrization/leakage 
scheme. In their work, they showed the low-|T/W| instability appears when the progenitor rotates 
2-7T rad/s or faster and alters the gravitational wave radiation. As for the explosion dynamics, 
they reported that bipolar outflow is driven by strong magnetic field if the initial central magnetic 
field strength is the order of > O(10 12 )G (|Mikami et alJ boQsl : IScheidegger et all boOfll ) ■ This is 
because the magnetic field of the order of > 0(1O 12 )G is easily amplified, a factor of ~ 0(1O 3-4 ) 
(jBurrows et al.l 120071 ). simply by the compression and the rotational winding effects during core 
collapse. Magnetic pressure with B ~ 10 15 ~ 16 G is comparable to matter pressure in the vicinity 
of the surface of PNS and thus drives bipolar outflow. The amplification factor (~ 0(1O 3 ~ 4 )) 
hardly depends on the initial strength unless some other nonlinear amplification mechanisms work. 
In the 3D context however the non axisymmetric motion may trigger the previously mentioned 
nonlinear amplification mechanisms (e.g., MRI) and may alter the amplification factor or the time 
scale comparing to axisymmetric motion. Therefore, especially when the initial magnetic field is 
much weaker than ~ 0(1O 12 )G, the axisymmetric assumption may not be suitable for the CCSNe. 
Additionally, since the strength of the order of ~ 0(1O 12 )G is considered to be unrealistically strong 
for pre-collapse stage, it should be examined how the initially weak (or more realistic) magnetic 
field, e.g., ~ 0(1O 9 )G, is amplified and whether it affects the explosion dynamics or not. 

Another challenge for numerical simulations of CCSNe is the treatment of general relativistic 
effects. Since the gravity plays intrinsic role and also some very massive progenitors (> 25M , 



see, iTanaka et al.l (12009a)) form black halls (BHs), we cannot say any conclusive results about the 
explosion mechanisms of CCSNe without taking account of the general relativity (GR). As same as 
the context of 3D MHD works of CCSNe, the r e are n ot so many works done by MHD simulations 
including GR effects (see, e.g., IShibata et al.l (|2006l ); ICerda-Duran et al.l (|2007l )) and furthermore 
3D GRMHD works of CCSNe have not been done yet. 

In this paper, we describe our newly developed two 3D-MHD codes for CCSNe simulation. 
Their features are as follows. Both codes employ adaptive mesh refinement (AMR) technique 
and can cover wide dynamical ranges from the central compact object (~ 10km) to beyond the 
several times of iron core radius (~ 10000km) or more. Self gravity is included in the Newtonian 
approximated MHD code and the dynamical metric is included in GRMHD code. High resolution 
shock capturing scheme is adopted in both codes and can handle the shock discontinuity without 
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numerical viscosity. We use the staggered mesh algorithm and adopt the Constrained- Transport 
method to evolve the magnetic field. We have done many test problems and confirmed their 
abilities. 

This paper is organized as follows. From Sec. [2]to[U we describe our numerical codes, their 
detailed properties and adopted techniques. In Secj5] and El we report several numerical tests of 
Newtonian approximated and general relativistic codes, respectively. Section [7] shows our first core- 
collapse supernovae calculations and their initial setups. We summarize in section [HJ We adopt 
cgs units for NMHD code and geometrical units for GRMHD code in which G = c = 1. In our 
practical calculations of a collapse of 15M star, in Sec. physical quantities are notated in cgs 
units. Greek/Latin indices run through 0-3/1-3. 



NMHD Code 



Our NMHD code solves ideal NMHD equations written in conservative forms together with 
source terms of self gravity, which are written in the following form: 

dQ 



at 



+ V-F = S 



and Poisson's equation for gravitational potential 



V 2 $ = 4vrGp 



(1) 



(2) 



Where, in Eq. 
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Eq. (H]) represents the conservation of mass, momentum, energy and the Faraday's law. p, u, B, 
E, Ptot, $ and I are rest-mass density, fluid velocity, magnetic field, energy density, total pressure, 
gravitational potential and 3x3 unit matrix, respectively. The energy density and the total pres- 
sure are expressed as 
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E = pe + 



+ 



2 8vr 
B 2 



tot 



p + 



S7T 



(6) 
(7) 



Here, p and e are gas pressure and internal energy density, respectively, and are related via the 
equation of state, p = p(p, e). We also define primitive variables P=(p, u, e, B) which is uniquely 
obtained from conservative variables Q via EOS. We employ staggered mesh and all variables 
except magnetic field are defined at cell center while, for instance, S x (i,j,k) is defined in cell surface 
(i+1/2, j, k). Hereafter we use B' = B/v4tt and do not express " ' ", unless otherwise stated. 

To solve time dependent equation (pQ), in our cartesian Eulerian grid, we ad opt Roe-type upwin d 
solver in which numerical flux P is defined in cell surface (see Fig(T]). Following iPowell et al.l ([1999), 
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Fig. 1. — Positions where numerical flux and electric field are defined. For instances, x directional 
numerical flux F x is defined at (i+1/2, j, k) and x component of electric field Ef x is defined at (i, 
j+1/2, k+1/2). 



F is constructed from the spectral decomposition of the system and expressed as 



(F(P L ) + F(P fl )) 



m=l 



(Qr - Qi)|A m |Rr, 



/2 



(8) 



Here the index L/R represents position immediate left /right of the cell boundary, i.e., when we 
evaluate the numerical flux at (i+1/2, j, k), the position L and R stand for (i+ 1/2-0, j, k) and 
(i+1/2+0, j, k), respectively. L m /R m and X m correspond to left/right ei genvectors and eigen values 
of the system, respectively, and defined in cell surface. Originally, in IPowell et al.l (|1999l ). they 
decompose the system into eight spectral modes in which one mode carries the monopole moment 
of magnetic field. Meanwhile, we adopt the constrained transport method for time evolution of 
magnetic field and we thu s consider only s even modes in our NMHD code (for explicit expressions 



of seven eigenvectors, see (|Rvu et al.lll995l )). Eigenvalues \ m are defined as 



\ m = (U — Cf, U-C_A, U-C s , U, U + C s , U + C4, U + Cf) 



(9) 
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Here u is the fluid velocity normal to cell boundary; cy the fast magnetosonic speed; c s the slow 
magnetosonic speed; cj± the Alfven speed, ca, Cf, c s are expressed by 



c f,s 



bJVp 



a 2 + B 2 /p ± J(a 2 + B 2 /p) 2 - 4a 2 c 2 4 



(10) 



(11) 



Where, f/s takes +/— in Eq. (jlip . U n represents magnetic field component normal to cell boundary 
and a is the speed of sound. The speed of sound a in a general form of EOS can be defined by 



P/P 1 



dp 



0e_ 

op 



(12) 



Since we define the conservative variables other than magnetic field at cell center, we have to 
interpolate to obtain those values at immediate lef t /right of cell bo undary. As for the interpolation, 
we adopt the monotonized central (MC) method (jVan Leerlll977l ) as 



i±l/2=F0 



MC(x,y) 



P i ±MC(P m -P i ,P i -P i _i) 

/ 2 sign(x) min(|x|, \y\, \x + y|/4) 




for xy > 
otherwise 



(13) 



Additionally, we empl oy total var i ation al diminishing (TVD) scheme which has second order con- 
vergence in space (see iRyu et al.l (|1995l )) and we briefly summarize how we implement it into our 
NMHD code below. In TVD scheme, summation respect to the spectral modes in Eq. ([8]) is 
modified as 



^ Lm • (Qr ~ Qi)|A m |R m = ^ (Ql (^~^ m + ^ m ) Ld Q™,i ~ (9l + 9r) J R-m 

m=l m=l V V x / / 

In Eq. (fT4|) we consider x directional numerical flux F x ,i+l/2 defined at (i+1/2, j, k) and 



(14) 
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At 
As 



sign(fifi+i) max[0, min(|g i+ i|, gisign(g i+1 ))} 

sign^i) max[0, mm(|#i|, ^-isign^i))] 

[ (9r ~ 9i)/LdQ m:i for LdQ m>i / 







otherwise 



+ e for \x\ < 2e 
\x\ for Ixl > 2e 



(15) 
(16) 

(17) 
(18) 

(19) 
(20) 
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/ 0.01 for m = 1,7 

0.1 for m = 2, 6 

0.1 for m = 3, 5 

V for m = 4 



(21) 



After we obtain the numerical fluxes in all sides of cell, the co nservative variab les other than 
magnetic field are updated through predictor and corrector steps (jPress et al.lll992l ). In predictor 
step, Q n is updated from time level n to n + 1/2 by At/2, i.e., 



Q 



71+1/2 



Q n +0.5Ai 
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pm pin 
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Az 




and in corrector step from time level n to n + 1 by At with using predicted values, i.e., 



Q 
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Q n + At 
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c n+l/2 



(23) 



This method is second order convergence with respect to time. 



On the other hand, as for the time evolut ion of magnetic field B, we adopt the constrained 
transport (CT) scheme (jBalsara &: Spicerlll999l ) in which B is evolved as 



<9B 

~dt 



V x E 



(24) 



Electric field E for CT scheme is defined at the cell edge (see Fig(T]) and is evaluated from Roe-type 
numerical flux, Eq. (151) , with appropriate interpolation. We simply express the electric field as 



E 



X,i,j+l/2,k+l/2 



pi6 _i_ pi6 

r y,i,j+l/2,k r y,i,j+l/2,k+l 



? 7 

z,i,j,k+l/2 



P 7 

z,i,j+l,k+l/2 



/4 (25) 



Here the upper suffixes in the right hand side denote components of the numerical flux, y and 
z components of the electric field are obtained straightforwardly by permutation as x — > y, y — > 
z, z -+ x and x — > z, y — > x, z — ^ y, respectively. 

In self gravitating system, the source term contains gravitational potential which is obtained by 
solving Poisson equation (our method to solve Poisson equation with AMR framework is described 
in Sec. 14. 2p . Since solving Pois son equation is very time consuming task, we solve it only once in 
one hydro dynamical time step (jTruelove et al.lll998l ) after predictor step is completed. Practically, 
in predictor step, we extrapolate gravitational potential <& n at nth time level via 3> n = (SQ^ 1 / 2 — 
<3? n_3 / 2 )/2 and then predict variables at (n + l/2)th time level through Eq. ()22|) . After the predictor 
step is completed, we solve Poisson equation by using p n+1 / 2 . In corrector step, we fully evolve 
the nth time level to (n + l)th variables by using Q n+1 / 2 and $ n+1 / 2 . Even though we extrapolate 
gravitational potential at nth time level, our numerical results do not show any large error in 
conservation of energy and numerical convergence is also achieved (described later in Sec. 15. 3ft . 
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On a final note, we mention about a numerical instability which is characteristic to Roe-type 
scheme. Even though Roe-type numerical flux is less numerically dissipative and has good shock 
capturing ability, one problem arises which is a so-called "odd-even decoupling". This instability 
appears when the shock normal is directed parallel to the gr id alignment. To avo id this instability, 
we adopt " carbuncle cure" in our NMHD code by following lHanawa et al.l (|2008l ) . 



GRMHD Code 



Formalism of our GRMHD code is based mainly on IShibata fc Sekiguchil (|2005l ) . It can be di- 
vided into two parts, one is MHD part and the other is Einstein's equation part. MHD part describes 
time evolution of matter on the background of spacetime metric and the metric is evolved accord- 
ing to Einstein's eq uation through the so-called "Baumgarte-Shapiro-Shibata-Nakamu ra (BSSN)" 



formalism (see, e.g.. IShibata &i Nakamura 



1995 



Baumgarte fc Shapirolll999l : lYo et all 12002^ . 



Before going to brief summary of our method, we describe our fundamental variables of MHD 
and metric parts. We set fundamental variables of MHD part as rest-mass density p, specific 
internal energy e, 4-velocity and magnetic field measured by a comoving observer. For metric 
part, the 3- metric 7^ and the extrinsic curvature Kij are th e fundamental ones a dopting the 3+1 
formulation of " Arnowitt-Deser-Misner (ADM)" formalism (jArnowitt et al.lll962l ). Then the line 
element of the spacetime can be expressed as 



ds 2 



-a 2 dt 2 + - fij (dx i + f5 i dt)(dx j + (5 j dt) 



(26) 



Here a, f3 l are the lapse function and the shift vector, respectively, and determined by arbitrary 
chosen gauge condition (Sec. I3.3() . Hypersurface of constant t is foliated in the spacetime so that 
a unit normal vector n^(n^) to this hypersurface becomes 



rr 



(1, -^)/a & n M = (-a,0) 



(27) 



Fundamental metrics 7^, are converted to 5 variables in BSSN formalism which are; the 
conformal exponent <j) = ln(7)/12, here 7 is the determinant of 3-metric 7^; the conformal 3- metric 
Jij = the trace of the extrinsic curvature K = tr (Kij); the tracefree extrinsic curvature 

Ay = e~ A ^{Kij — yijK/3); and three auxiliary variables F$ = 5 3 Tijfc, here ",k" represents partial 
derivative with respect to k direction. Hereafter, Di and Di denote the covariant derivatives with 
respect to 7^ and 7^, respectively. 

Stress energy tensor for ideal magneto-hydrodynamical fluid is expressed as 

T^u = {ph + b 2 )u^u v + {p + b 2 /2)g^ - b^b v (28) 

Here h(= 1 + e + p/p) is enthalpy and we define magnetic pressure as b 2 /2 = b^b^/2. With these 
settings, we define other useful quantities which are; 3-velocity v 1 observed by Eulerian observer at 



rest; and magnetic field observed in the fluid flame, i.e., = n^. 

v l = u l /u l 

B» = (0,e 6< ^(W-a&V)) 
Here W = au l is the Lorentz factor. We also define primitive/conservative variables P/Q as 



(p, Ui,e,B l 



Q 



p* 
Si 




T 




B l 





pe^W 
e^((ph + b 2 )Wui- ab%) 
e^iiph + b 2 )W 2 - (p + b 2 /2) - {ab 1 ) 2 ) - p % 
B l 



(29) 
(30) 

(31) 
(32) 



3.1. Magneto Hydro dynamical Equations 



Basic equations of magneto-hydrodynamical part in general relativistic form are written as the 
following conservative- like equations. 



dtp* + di(p*v % ) = 



(33) 



d t Si + O^Sivi + ae 6<t> P tot 5j - B^h/u 1 ) = 
-S d ia + S k d^ k + 2ae^S k k 8 l <t ) - ae 2 *(S jk - P to t7ifc /2 



(34) 



d t r + diiSov* + e 6 ^P tot (v l + (3 l ) - ab'B^u 1 



p*v 



(35) 



d t B i + dj (BV -v i B j ) = (36) 
Where So = r + p* and these equations can be expressed in the form of dtQ + di~F % = S. 

Our procedure to evolve the MHD conservative variables Q i s sim ilar to our NMHD code 
except we use the HLL (Harten-Lax-van Leer) flux (jHarten et al.l Il983l ) and not Roe-type one. 
HLL-flux is less numerically expensive compared to Roe-flux, since we only have to consider the 
two fastest left and right going wave speed without considering eigenvectors like expressed in Eq. 
([8]), however HLL-flux has sufficient capability to follow shocks and is suitable for our aims. The 
fastest left /righ t going wave speed are evaluated from 4th order eigenvalue problem (see, Eqs. (58), 
and (63-65), in I Anton et al.1 120061') . We solve this p roblem by iterative Newton method with given 
sound speed. Following IShibata h Sekiguchil (|2005h . the sound speed a is defined as 
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By solving 4th order eigenvalue problem the fastest left/right going wave speed A1/A7 (correspond 
to A1/A7 in Eq. ([9])) can be obtained. We evaluate Ai and A7 at immediate left and right of the cell 
boundary by adopting several types of reconstruction schemes such as monotonized central (MC) 
or piecewise linear method (PLM). In this paper, we adopt only MC method which has second 
order convergence with respect to space. Then the fastest left/right going wave speed A_/A + at 
cell boundary are defined as 

A_ = max(0, Ai )L , \ 1>R ) (38) 
A+ = max(0, A 7>L , A 7)jR ) (39) 



With these wave speed, we define HLL flux, by following lAnton et al.1 (120061 ). as 



F A + F(P L ) - A_F(P fl ) + A„A + (Q fi - Ql) , . 

* HLL = ~ = (40) 

A+ — A_ 

Here A = A/a and F(P) is an appropriate flux vector in Eqs. (|33ti36|) . Solenoidal constraint of 
magnetic field is satisfied by CT scheme as the same procedure as NMHD. 

Once we update conservative variables Q, we have to obtain primitive variables P by solving 
following three coupled equations with iterative Newton method. 

r = t( P ,s,W) (41) 
S l S t = S 2 (p,e,W) (42) 
p* = p*(p,W) (43) 

We employ the same recovering procedure as proposed in Cerda-Duran et al. (|2008l ) with adopting 
"safe-guess values" when the iteration does not converge. 



3.2. The BSSN Equations 



Next we describe our method to evolve metric part. As previously menti oned we evolve BSSN 
variab l es (</>, 7,;,-, K, Ajj, FA acco r ding to followin g equations (see, e.g., IShibata h Nakamura 
(|l99.^ : lBaumgarte fc Shapiro! d 19991 ): I Yo et al.l feood )). . 



(d t - Cp)Aij 
(8 t - C P )K 

(d t - p k d k )Fi 



-2aAij 
-aK/6 



DiDjaf* + a(KA tJ - 2A ik ^ kl A 



Aa + u(AijA ij + K 2 /3) + 4vra(,Soe" 6 * + f j S i:j ) 



= -ISirae-^ Si 

+ 2a [f jk A ij>k + ffrfiAij - & k h jk>i /2 + 6<f> tj A{ - 2K A /3 
+ 5 jk \-2a, k Aij + f3\ k h lU + (i u P l j + ~ 2j lJ p\ l /3) tk 



(44) 
(45) 
(46) 
(47) 



(48) 
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T 



8]), Cp is Lie derivative with respect to f3 l ; "trf" denotes trace-free operator; 
S lJ and hij = 7^ — bij . Rij is the Ricci tensor and consisted of two parts in 



In these Eqs. 
A = D i D i ; P 
the form of 

Rij = Rij + Rf- (49) 

For explicit forms of the Ricci tensor and several notes when calculating the Ricci scalar, see 
Shibata &: Urvul (120021) . W e evolve these BSSN variables by the second order scheme in space (e.g., 



Appendix of IShibatal 1 19991 ) and by the iterative Crank-Nicholson scheme with three steps in time. 
M any recent numer i cal si mul ations in full general relativity adopt fourth order scheme in space such 



Zlochower et al.l (|2005l ) or lEtienne et al.l (|2008l ). However, such higher order scheme is necessary 



as 

especially when the metric is highly distorted such as around the BH. Our numerical simulations 
with a relatively low mass star (15M ) do not show any BH formation and we thus consider second 
order scheme is acceptable at this time. 

For the metric, there are several mathematical and physical constraints. As for the mathemat- 
ical constraints, det(7ij) = 1 and tr (Aij) = should be satisfied. We enforce following two artificial 
procedures 



lij 
A 



'.1 



7ij/det(7y) 

Aij - j kl A k aij/3 



(50) 
(51) 



after each update to maintain numerical stability. Physical constraints are the Hamiltonian and 
momentum constraints. 



n 



Mi 



~ ■ ~ , e^R e 5 ^ 
D l D^ - i-p + 2vrS e-^ + — 



Aij A - 1 



-K' 







Dj ( e^Al, 



-e^DiK - 8nSi = 



(52) 
(53) 



We do not enforce any artificial modifications to satisfy these constraints, though monitor these 
values just as our code check. However, we enforce Hamiltonian constraint every time we re- 
fine/coarsen the AMR blocks by solving above Poisson like non-linear equation(|52|). We monitor 
Cu defined by 

1 f pM 



\D i D i e'l > \ + 



\e<t>R\ 



+ |27r5 e-^| + 



■dx 6 



(54) 



to check the Hamiltonian constraint and accuracy of our code. Here, Mb ar is proper rest mass and 
defined by Eq. ([5TD . 



3.3. Gauge Conditions 

Hyper-surface at constant time t can be foliated in spacetime arbitrary, but is usually chosen 
so as to keep time evolution numerically most stable. As for the time slicing condition which deter- 
mines lapse (a), we adopt several choices in our GRMHD code such as "the approximate maximal 
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slicing" or "harmonic slicing" or "1+log" slicing conditions. Since the approximate maximal slicing 
condition requires to solve poisson like non-linear equation every time step and very time consuming 
method, we usually adopt 1+log gauge condition given by 



d ta = ftdia - 2aK 



(55) 



We have implemented "dynamical gauge condition" for shift (/3*), following IShibatal (|2003l ). which 
is given by solving 

d t (3 l = fi (Fj + AtdtFj) (56) 

Where, At is the numerical time step. By imposing these gauge conditions, we do not suffer from 
time consuming Poisson like equations at every time step, while we do not encounter any numerical 
instabilities throughout our calculations of core collapse of massive star. 



3.4. Diagnostics in GRMHD 



In GRMHD code, global quantities such as total baryon rest mass Mb ar ; ADM mass Madm; 
total angular momentum along the rotational (z) axis J z ; internal energy -Ejnt j magnetic energy 
E mas ; kinetic energy T^; rotationa l kinetic energy T rot are defined as expressed below (see, e.g., 
Duez et alJliool iKiuchi et aJbood ). 



M bar 
-Madm 
J, 

^kin 
Trot 



p*dx 3 



=5</> 



S ° + lfor" \ AijA%J ~ 3 " f3Rij€ 



-ic6 



dx" 



mag 



AX _ AV 

— y — h xS y - yS x + 

p*edx 3 
pJiUiV % dx 3 

pjiu^dx 3 

b 2 
2 



xK, y - yK^ x xj^Aij - yj'xAij 



12vr 



16vr 



— We^dx 3 



(57) 
(58) 

(59) 

(60) 
(61) 
(62) 
(63) 



Then the gravitational potential energy E^ v is defined by Egj- V = — (Mb ar — Madm + E{ nt + 
T]j m + E mag ). Because of our formulae for MHD part and because of our AMR scheme (see, Sec. 
I4.3p . Mb ar is conserved with high accuracy. On the other hand, conservation of Madm which is 
guaranteed from the Einstein's equation in the absence of gravitational radiation is violated due 
to the accumulation of numerical errors in our CCSNe simulations. Then, several % fluctuation 
appears (see, Fig{16]in Sec. [6]) which is approximately 2-3 orders larger compared to that of Mb ar 
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in our CCSNe simulations. Consequently even though our initial condition satisfies E grv < 
(i.e., gravitationally trapped system), there sometimes appear that E gTV > during calculation. 
Therefore we estimate the gravitational potential energy as E gTV ~ —(-Bint + ^kin + -^mag) and apply 
it to such as the rotational to gravitational energy /3 rot = T ro t/|^ g rv|- 



4. Adaptive Mesh Refinement 



One of difficulties in computational astrophysics is that we have to handle wide dynamical 
range in a limited computational resource. For instance, in the context of CCSN simulation, the 
PNS is a size of ~10 km and on the other hand radius of the iron core is the order of > O(1000)km. 
If we cover such a vast range (several times of the iron core) with a uniform resolution, e.g., ~ 500m 
to resolve interior of the PNS, it becomes impossible to calculate with our limited computational 
resource. We thus raise resolution in the vicinity of proto-neutronstar and, at the same time, lower 
resolution far from the centre. To realize such situation, we incorporate the AMR technique (e.g., 



Berger Colella 



1989 



Powell et al.lll999l ) into our codes. 



4.1. AMR Structure 

In our codes, computational domain is divided into "blocks" (hereafter, AMR block) and every 
AMR block consists of 8 x 8 x 8 cubic cells and of 2 additional cells as ghost zones in every side 
of block. Every AMR block belongs to a refinement level "I" and if the refinement level I is raised 
by one, the AMR block is divided into 8 blocks with halved cell width. On the other hand, if all 
2x2x2 neighboring AMR blocks are assigned to lower their refinement level by one, 8 blocks 
merge into one block with twice cell width. 

Our codes are fully parallelized and adopt Message Passing Interface (MPI) for communication 
between different nodes. Then it requires load balancing in AMR frame work. Our method for this 
purpose is like this. Three dimensional structure of AMR blocks are projected on one dimensional 



structure connected by "Hilberf space filling curve (jHilb ertl 1 1 89 ll ) . Along the Hilbert curve, AMR 
blocks are numbered sequentially. Then AMR blocks projected on one dimension are allocated to all 
computational nodes in a straight forward manner. Using such a projection scheme, e.g., connecting 
by Hilbert curve, enables us to minimize the data transfer between different computational nodes. 
This is because surface area of a "chunk" of AMR blocks allocated in one node by above method is 
minimized as much as possible and, thus, we can minimize time to spare for the data communication. 
In top-left panel of FigJH we display an example of Hilbert curve in two dimension. In this 
figure, asterisks(circles, crosses) denote centers of cells with width 0.5(0.25, 0.125) and lines are 
Hilbert curves. As seen in this panel, two dimensional structure of AMR blocks is projected on 
one dimension. In right-top panel, we again display Hilbert curve which fills two dimensional 
computational domain covered by blocks of four different AMR levels. Background colors represent 
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Fig. 2. — Top-left: Hilbert curves filling two dimensional area covered by uniform AMR blocks 
with three different levels. Top-right: Hilbert curves filling two dimensional area covered by 
nonuniform AMR blocks with four different AMR levels. In this panel, AMR blocks are allocated 
to 8 nodes and blocks allocated to one node are connected through one sequential line. Bottom: 
Three dimensional extension of top-left panel with two different levels. 

AMR levels and blocks connected one curve are allocated to one computational node. Therefore, 
in this panel, all AMR blocks are allocated to 8 nodes with maintaining load balancing. In bottom 
panel, we also display three dimensional extension of Hilbert curves with two different AMR levels 
for reference. 

By adopting such method, our AMR structure has flexibility to refine or coarsen AMR blocks 
locally. 



4.2. Poisson Solver under the AMR Framework 

In NMHD code, we have to solve Poisson equation in the form of Ax = B for the self gravity 
and; in GRMHD code, Poisson like non linear equation in the form of Ax = B(x) for the initial 
Hamiltonian and momentum constraints. Here A is a given (A^iock x 8 x 8 x 8) x (Abiock x 8 x 
8x8) matrix, B is a given (Abi oc k x 8 x 8 x 8) vector and x is a solution we seek which is 
(iVbiock x 8 x 8 x 8) vector. Abiock is a total number of AMR blocks. In GRMHD, B(x) is a 
vector containing non-linear term o f x. We adopt itera ting method, the so called "BiConjugate 



Gradient Stabilized (BiCGSTAB)" (jvan der Vorstlll992l ) method to solve such huge simultaneous 



equations. Our strategy for solving this equation under our AMR structure is; (1) we set an AMR 
level I which is at initial; (2) for all AMR boxes whose AMR levels are larger or equal to I, we 
project their physical quantities, such as the density, to boxes of AMR level I by the coarsening 
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procedure and construct A Sz B; (3) we then solve equation by BiCGSTAB on the uniform mesh 
with appropriate boundary conditions at the interface of AMR level I and I — 1 and also at the 
outer boundary; (4) increment I by one and repeat these procedures from (1) again. Note that, 
we have to pay special attention at the interface of different AMR level, i.e., I and I — 1. At 
here, we have to connect both the solutions and their first derivation smoothly, otherwise there 
appear some non- physical divergence . To avoid this, we adopt quadratic and bilinear interpolation 
methods following iMatsumoto I (|2007l ) to evaluate ghost zone values of AMR level I, e.g., <&(i — 1, j) 
and &(i — 2, j) in FigJ3J Here, we summarize our interpolation method in two dimension. Three 
dimensional extension can be done in a straightforward manner. If we seek gravitational potential 
— 1, j), we first have to obtain &(A) at (I, J — 1/4) which is derived via 



*(A) = $(/, J) - ^MC(<&(/, J + 1) - $(/, J), $(J, J) - <&(/, J - 1)) 



(64) 



Here, MC(x,y) is the monotonized central method. Then — 1, j) is derived via quadratic 



J 



y 



i + 1 



x 



Fig. 3. — Schematic figure of interpolation of ghost zone value. Four small boxes and one large box 
belong to AMR level I and I — 1, respectively. 

interpolation 



10$(i,j) - 8$(A) + 3#(j + 1, j) 
15 



(65) 



<&(B) which is required to evaluate $(z — 2,j) is obtained by following equation 



= $(/,J) + $(J-l,J) _ 1 $(J,J + 1) + $(I-1,J + 1) _ $(J,J) + ^(/-l,J) 
1 j 2 2 1 2 2 



$(/, J) + *(/ - 1, J) $(/, J - 1) + $(/ - 1, J - 1] 



) (66) 



Three dimensional extension of this smoothening method is done by replacing, e.g, Eq, (|64p . with 
bilinear interpolation 



*(A) = $(I,J,K) - ~MC($(I,J+1,K)-<S>(I,J,K),$(I,J,K)-$(I,J-1,K)) 

- ^MC($(I, J,K + 1)-Q(I,J,K), J, if) -*(I,J,K-1)) (67) 
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On the other hand, lower level ghost zone value + 1, J) is derived by "restriction" procedure 
and this is simply averaging $ over the 2x2 (or 2 x 2 x 2 in 3D) adjacent cells via 

$(I+1,J) = - $(i + n,j + m) (68) 

n,m=0,l 



4.3. Boundary of AMR Blocks 

To guarantee the conservation law and the solenoidal constraint of magnetic field, we have to 
reflux the numerical flux, Eqs. ([8]) and (|40p . and the electric field, Eq. (|25p . at where AMR boxes 
of different levels are contacting. In FigJH we display schematic picture of refluxing procedures. 
For instance, the numerical flux F x belonging to AMR level I — 1 and defined at cell boundary is 
replaced by summation of f x ,i(i = 1,2,3,4) which belong to AMR level I. Similarly, as for the 
electric filed, the electric field E y defined at cell edge is replaced by summation of e y ^{i = 1,2). 
These procedures ensure the conservation and the solenoidal constraint below the round off error. 




Fig. 4. — Schematic picture of refluxing of the numerical flux F x (f x ) and the electric field E y {e y ). 
Left two and right boxes belong to AMR level I — 1 and I, respectively. 

Additionally, as for the ghost zones, we have to obtain physical variables every after the time 
updating and this procedure is sorted into three cases. 

(1) For AMR box whose neighbor has the same AMR level, we simply copy all physical variables. 

(2) For AMR box whose neighbor has lower level, ghost zone variables are interpolated and are 
sent from lower to higher level box. For this interpolation, we use the same method as used in our 
Poisson solver (described in Sec. 14. 2D other than the magnetic field B. As for the magnetic field, 
we have to interpolate whi le maintaining the solenoidal constraint and this is done by adopting the 



same method proposed by iBalsaral (|200ll ). 



(3)If AMR level of the neighbor is higher, physical variables are evaluated by "restriction" 
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procedure. In our codes, this restriction procedure is simply averaging the variables of 2 x 2 x 2 
adjacent cells which is the same method as used in our Poisson solver (Sec. I4.2|) . 



4.4. The BSSN Evolution Under the AMR Framework 

During time evolution of the BSSN variables, we have to derive the spatial derivatives of 
metrics not only along one direction, e.g., d x , but also the cross derivatives, e.g., d xy to obtain such 
as the Ricci tensor. Therefore, if there exist some discontinuity in the spatial derivatives across the 
AMR refinement boundary, spurious oscillations of the BSSN variables appear near the refinement 
boundary. Since the time marching is simultaneous across all the AMR boxes in our codes, there 
is no time lag between different AMR levels. However, we should carefully interpolate the buffer 
zone's metrics especially for AMR boxes whose neighbors have lower AMR levels than theirs. This 
situation is the same as that appeared in our Poisson solver (Sec. 14. 2p and we adopt the same 
strategy to obtain the buffer zone's variables. In addition we have to evaluate the metrics along 
the edge of AMR block (e.g., (i,j,k) = (9,9,1 ~ 8)) for the cross derivatives and this evaluation 
can be done in a similar manner as that in the normal buffer zone's case. For instance, if we seek 
a metric X at (i,j,k) = (9,9, 1), X(A) which corresponds to $(A) in Eq. (j64l) is replaced by 

X(A) = X(l, 1, 1) - ^MC(X(1, 1, 2) - X(l, 1, 1), X(l, 1, 1) - X(l, 1, 0)) (69) 

Then, in a straight forward manner of Eq. (|65p. X(9, 9, 1) can be derived by 

X(9, 9, 1) = 1°A-(8,M)-8*(A) + 3 A-(7.7,1) (fl)) 

However, we cannot completely suppress the spurious oscillation s in the refinement b oundary 



and, in such case, adding numerical dissipation is sometimes useful (jSchnetter et al.ll2004l ). Even 
though we do not add the dissipation at present codes, we do not suffer from growth of the noises 
and the code crash. Several tests of the BSSN evolution with AMR structure are summarized in 
SecEl 



5. Tests for NMHD Code 

In this section, we introduce several test problems done by our newly developed 3DMHD code 
in the Newtonian approximation. We verified the accuracies of our new AMR-NMHD code through 
several test suites. 
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5.1. One Dimensional Shock Tube Test 



We calculated several 1-D hydrodynamical shock tube tests and com pared to the exact solutions 
obtained from the code HE-E1RPEXACT of the library NUMERIC A (jjorcj \l9M) . We show one 
test in Figj5]in which we assumed ideal gas with adiabatic index 7 = 1.4. The mesh width of exact 
solution is 1/1000, meanwhile, the "effective" mesh width of our numerical run's are 1/2048, 1/512, 
1/128 for blue, green and red dots, respectively, 
test 1) 



tot J 



(1.0, 0.0, 1000.0) x < 0.5 
(1.0, 0.0, 0.01) x > 0.5 



For MHD shock tube test, we show the same test described in (|Brio fc Wul [1988) with three 





Fig. 5. — Results of test 1) at t = 0.012. Black solid line is analytical solution and blue, green 
and red dots/lines are our numerical results. Each color represents the maximum refinement level 
-^AMR.max- We set LAMR.max = 2, 4, 6 for red, green and blue, respectively. 

different types of the EOSs in FigEJ 
test 2) 

(p , V , B x , B y , B z , P tot ) = 

f (1.0, 0, 0.75^4^, 
| (0.125, 0, 0.75\/47, 

Initial total pressure -Ptot) which excludes the magnetic pressure, is divided into three parts such as 
the gas (-Pgas)j the degenerate (Pdege) and the radiation (P ra d) pressure term. 



4vr, 0.0, 1.0) x < 0.5 
4tt, 0.0, 0.1) x > 0.5 



a) : Ptot = Pgas = pT 

b) : Ptot = Pgas + Pdcge = pT + 0.3p 4/3 

C) : Ptot = Pgas + Pdcge + Prad = pT + 0.3p 4 / 3 + T 4 /3 

Internal energy has also three parts corresponding to the gas (e sas ), the Fermi (ffdege) and the 
radiation (e ra d) energy as defined by the following. 

£tot {Pi T) = £ g as + £dege + ^rad 
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tin-, Pdcwa , , 3-P ra( j 

+ / — 7^dp-\ 



(7 - 1)P Jo P 



Here 7 = 2.0 and EOS a) corresponds to the original model reported in lBrio fe Wul ()1988l ). From 



P=pT+0.3p 4/3 +T 4 /3 
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Fig. 6. — Results of test 2) at t = 0.1. Each line corresponds to different EOS shown at the top. 



FiglMSl we see that our code captures the discontinuities accurately and also no numerical insta- 
bilities are seen. 



5.2. Poisson Solver 

As for the tests of our Poisson solver, we set two types of spherically symmetric density 
distribution like below which have analytical solutions and compare our results with analytical 
ones. 

test 3) Homogeneous sphere of radius R and density po 

test 4) Centrally condensed sphere with density distribution p(r) 



po 



r < R 



P(r) = I (71) 

r > R 



These tests are the same tests done in lStone &: Normanl ([1992]) and have analyti cal solutions, thus 



we ca n easily check the accuracy of our Poisson solver (for analytical formulae, see lStone &: Norman 



19921 ). We set R = 10 8 cm, p = 10 12 g cm" 3 and r c = 3 x 10 7 cm for both tests 3) & 4). Figd 
shows our numerical results with comparing analytical ones. Upper and lower two panels show 
|dV0| = |(V0ana - V0)/V0 ana | and |d</>| = |((/> ana - 0)/0 ana |, respectively. Here, <j) is the numerical 
result of the Poisson equation and ana is the analytical one. Deviations of our numerical results 
from the analytical ones are ~ 0.1% and we also find neither kink nor jump of both <j> and V</> at 
the interface of different AMR level boxes. Therefore, we consider our Poisson solver under AMR 
structure works with sufficient accuracy. 
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Homogeneous Sphere 



Condensed Sphere 




Fig. 7. — Results of test 3) (left) and 4) (right). See text for the definitions of 
Five different levels of AMR boxes are drawn. 



and |dV0|. 



5.3. Energy and Angular Momentum Conservations 

In this subsection, we check our NMHD code's accuracy against the energy and the angular 
momentum conservations. Since we consider one of energy source of the formation of bipolar outflow 
is the extracted angular momentum, we have to carefully trace the time evolution of angular 
momentum. As for the test of angular momentum transfer, we follow the collapse of a non- 
magnetized and rotating 25M star with the adiabatic gas with index 7 = 1.4. If no magnetic 
field exists and the fluid is adiabatic gas, the angular momentum is not transported. Result is 
shown in FigJHJ Abscissa and vertical axes represent the specific angular momentum and the 
total mass in solar mass unit which is the summation of fluid elements having less or equal to 
the corresponding specific angular momentum on the abscissa axis, respectively. If the angular 
momentum conservation is maintained, the curves do not change its form in time and we thus 
see the angular momentum conservation is well maintained from this figure. As for the energy 
conservation test, we calculate collapse of a rotating and magnetized 15M star (corresponds to 
model "NB12R020Sf" described later in Seed. I n Fig® we display time evolutions of various 
energy components and the error in energy conservation in left panel and the magnified view 
around the time of core bounce with different numerical resolutions to see the numerical convergence 
in right panel. As for the numerical convergence test, the computational domain is chosen as 
(x, y, z) = [—5000, 5000]km for solid curves and (x, y, z) = [—4000, 4000]km for dashed ones. Then 
the minimum grid widths become Ax m j n = 600(solid) and 480(dashed) m. From Fig® we see the 
energy conservation is well maintained and also see the numerical convergence is achieved within 
the range of our adopted numerical resolutions. 
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L, (cm 2 ."'J 



Fig. 8. — The angular momentum conservation. Abscissa and vertical axes represent the specific 
angular momentum and the total mass in solar mass unit which is the summation of fluid ele- 
ments having less or equal to the corresponding specific angular momentum on the abscissa axis, 
respectively. In the inner mini-panel, the time evolution of the central density is displayed. Each 
Color of line represents the elapsed time from the initiation of collapse and corresponds to the same 
color-coded circle in the mini-panel. 




50 100 150 60 65 70 75 80 85 

time (ms) time (ms) 



Fig. 9. — Left: The error in energy conservation and time evolutions of various energy components. 
E grv , Ei n t, Ekin and E mag represent gravitational, internal, kinetic and magnetic energy, respectively. 
Relative error is defined by a deviation of (E grv + Ej nt + Ek in + E mag )/|E grv | from its initial value 
and is shown in upper part. Energy conservation is maintained within ~0.2% error through our time 
evolution. Right: Magnified view around the time of core bounce to see the numerical convergence 
with respect to the grid resolution. Solid and dashed curves correspond to models with minimum 
resolution Ai m j n = 600m and Ax m ; n = 480m, respectively. The internal, kinetic and magnetic 
energies are normalized by 10 52 , 10 51 and 10 48 , respectively. 
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6. Tests for GRMHD Code 

In this section, we introduce several test problems done by our newly developed 3DGRMHD 
code. 



6.1. One Dimensional Shock Tube Test 



As for the basic test, we calc ulated the relativistic Brio & Wu MHD shock tube test (see, 
Brio fc Wulll988l ; lKomissarovlll999l ) with fixed flat metric and the results are shown in FigllOi In 




-0.4-0.2 0.2 0.4 

position 



-0.4-0.2 0.2 0.4 

position 



Exact 

Numerical 

£«,/? )=(2,0.4) + 



Fig. 10. — T he relativistic version of Brio &; Wu MHD shock tube test. Exact solutions are 



obtained from iGiacomazzo &: Rezzollal (|2006l ) and drawn with blue lines. T is the Lorentz factor 
and other notations are the same as in Figj6j Numerical results (red lines) are models with the 
gauge condition (a, j3 x ) = (1.0, 0.0) and black crosses are the one with (a, f3 x ) = (2.0, 0.4). Levels 
of AMR blocks are from to 3 and the highest resolution is Ax = 1/512. Time slice is chosen at 
t = 0.25 for a = 1 models and t = 0.125 for a = 2 model. Results are shifted to x — > x + 0.4 x t 
for non-zero shift gauge model. 



this figure, we show two models with dif f erent g auge conditions and compare with the exact solutions 
obtained from IGiacomazzo Rezzollal (|2006l ). One is (a, (3 X ) = (1.0,0.0) and represented by red 
lines and the other is (a, f3 x ) = (2.0,0.4) and plotted by crosses. Time slice is taken at t = 0.25 
and t = 0.125 for a = 1 and 2, respectively. In addition, for non-zero shift gauge ((3 X = 0.4) model, 
results are shifted to x — > x + 0.4 x t to let it coincides with other results. We can see our GRMHD 
code can handle the non-zero gauge conditions and shocks. 



6.2. Bondi Accretion 



In this subsection, we test our GRMHD code in a strongly curved and fixed spacetime. As for 
the test, we evolve the Bondi accretion flow with /with out magnetic fi eld a nd compare our results 
to analytical solution which are obtained according to lHawlev et al.l (I1984T). It is know n that the 
radial magnetic field does not influence the Bondi accretion ()De Villiers &: Hawlevll2003l ) and, thus, 
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we add the initial magnetic field via 



B = V x 



B n 



r(r + z) 



-y, x, 0) 



In this test, we adopt Kerr-Schild coordinate 



a 



In 



2 



2 + r J r 

1 + 2/r, r 2 , r 2 sin 2 6>) 



(72) 



(73) 

(74) 
(75) 



Here, the metric 'ju is written in spherical polar coordinate. Event horizon locates r = 2 and outer 
boundary is set at |a;,y, z| = 10. We excise the computational domain \x,y,z\ < 1.5 and simply 
connect the non- and excised region with first order extrapolation. We run four models with two 



different resolutions (N\ 



block 



8 3 , 16 3 ) and B 



r 

Horizon 



1 (or b 2 /p\ 



Horizon 



2.446) and ^ orizon = 0. 

In FigQTJ we show the rest density profiles of magnetized models at t = 100M in left panel and 
L2 norm of the errors in density in right panel. The errors in iVbiock = 8 3 models are multiplied by 
(1/2) 2 to show the second order convergence. From this test, we see that our GRMHD code which 
employs HLL flux with MC limiter actually reproduces the second order convergence and also that 
our code can treat the strongly curved spacetime. 





Fig. 11. — Left: Density plots of magnetized Bondi accretion test with different resolutions taken 
at t = 100M. Solid line represents analytical solution and the vertical dash-dotted line is where 
we excise the computational domain. Asterisks and crosses represent our results with different 
numerical resolutions. They correspond to iVbiock = 16 3 (Asterisk) and A^iock = 8 3 (cross), re- 
spectively. Right: L2 norms of the errors in density from analytical value. Long-dashed lines 
are non magnetized accretion models and solid lines arc magnetized, models with. B T \ Horizon — 

1.0 

(or 6 2 /p|Horizon = 2.446). Thick lines are iVbiock = 8 3 models, while thin lines are iVbiock = 16 3 
models. Note that the errors of -/Vbiock 
order convergence. 



8 3 models are multiplied by (1/2) 2 to show the second 
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6.3. div(B) = Constraint 

To check whether our two codes satisfy the solenoidal constraint for the magnetic field, we have 
checked the value div(B)/B. FigfTJ] display div(B)/B in the central [0, 100] x [0, 100] x [0, 100]km 3 
region from one representative model of GRMHD models. In this region, three different AMR 
level boxes exist. At this moment, it almost reaches core bounce time and three times refinement 
procedures have been done since the initiation of calculation. From this figure, we see that solenoidal 
constraint is well maintained below the round-off error which tells us that both the refinement 
procedures and the electric field refluxing work well. 

xicr' 4 



50 100 150 

R (km) 

Fig. 12.— div(B)/B is plotted of model "GB12R020Sf ' . In this figure, three different AMR levels 
are plotted in the central [0, 100] x [0, 100] x [0, 100]km 3 region. 




6.4. Test Problems With Dynamical Background 

In this subsection, we test our dynamical metric solver written by the BSSN formalism with 
and without matters. 



6.4-1- Linearized Teukolsky Wave 



First test is to follow the linearized gravi tational waves, the s o calle d "Teukolsky wave" 
(|Teukolskvl 1 1982 ) . in a vacuum space. Following IShibata fc Nakamural (|1995l ). we adopt the same 
mode I = \m\ = 2 and the initial wave amplitude is set to C = 0.01. For time slicing gauge 
condition, we adopt "1 + log" condition (see, Sec. I3.3p . In this test, we also check the influence 
of the boundary where different AMR level boxes are contacting. In FigJ13|, the initial condition 
(a component of the extrinsic curvature, K12) and AMR structure are displayed. As seen in this 
figure, the inner region of \x, y, z\ < 5 is covered by maximum AMR level meshes which we vary as 
Lamr = 1)2 and 3 to see convergence of the error. We extract the metric variables 7 TO and 7 2Z at 
(x, y, z)=(4.2, 0, 0) and compared them wi th analytical ones in FiefHl Red curves are analytical 



solutions (see, e.g., iNakamura et al.l (jl9871 )) and black curves are numerical results. Numerical 
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Fig. 13. — Our AMR structure and initial value of the extrinsic curvature Kyi in the equatorial 
plane of linearized "Teukolsky wave test". The inner region of \x,y, z\ < 5 is covered by higher 
resolution meshes. 

resolutions are Ax min = 0.078(solid), 0.156(dash-dotted) and 0.312(dashed). From top panel of 
Fig ll4l we see that the errors decrease with increasing numerical resolution. 
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Fig. 14. — Results of linearized "Teukolsky wave test" . Bottom and middle panels represent j zz — 1 
and 7 ra — 1, extracted at (x, y, z)=(4.2, 0, 0), respectively. Red lines are analytical solutions and 
solid, dash-dotted and dashed lines are results of maximum AMR level 3, 2 and 1, respectively. 
In model with maximum AMR level 3, the minimum grid resolution is Az m ; n = 0.078. Top panel 
represents deviations of ^ yy from analytical value with three different grid resolutions. 



6.4-2. Rotating Neutron Star 

Next test is an evolution of rigidly rotating neutronstar in equilibrium state. Initial parameters 
are the central rest density p c = 10~ 3 and the central angular velocity £l c = 10~ 2 . We use the 
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polytropic EOS P = Kff where K = 10 and T = 5/3 and assume rigid rotation. With these 
parameters, the central lapse is a = 0.701 and the ADM/baryon masses are Madm = l-49/Mt, ar = 
1.55. Outer boundary is taken at \x, y, z\ = L out = 17 and the equatorial radius of the NS is ~ 13.1. 
In FigJ15t we display 4 models with two different numerical resolutions Ax m i n = L out /32 (models 
A and C), Ai m ; n = L out /64 (models B and D). Models C and D are evolved with fixed matter 
distribution and only metrics are evolved, while, in models A and B, both matters and metrics 
are evolved. Upper and lower panels display time evolutions of deviations of the central lapse 
and the ADM mass from their initial values, respectively. In the lower panel, we also display the 
time evolutions of baryon mass for models A and B with thick lines. From this figure, we see that 
models C and D keep their initial configurations within 1%, while fully evolved models A and B 
show gradual decrease(increase) of ADM mass(central lapse). In this test, our treatment of the low 
density region outside of the NS is like this. We set the floor density value as Pfl 00 r = 10~ 9 x p c ,max 
and, in every time step, for all cells whose density is smaller than 10 x pq oot , we assume them 
as vacuum. Then their density and velocity are reset to /9fl 00 r and 0, respectively. We consider 
that this treatment is too simple and, e.g., the conservation of baryon mass is violated as seen in 
thick lines. However, in the context of CCSNe, we currently do not have to treat vacuum space 
and, additionally, models A and B show numerical convergence with respect to grid resolution, we 
consider our GR code works with sufficient level for our aims. 




Fig. 15. — Time evolutions of the central lapse (top) and the ADM mass (bottom) normalized by 
their initial values. A (dash-dotted) and B (dotted) are fully (both matters and metrics) evolved 
models, while, in C (dash) and D (solid), only metrics are evolved. In the bottom panel, we also 
plot the time evolutions of baryon mass with thick lines. Grid numbers are 64 3 for A and C and 
128 3 for B and D. 



6.4-3. Box Refinement and Numerical Convergence 

In the end, we mention about the influence of AMR refinement procedures during the collapse 
and about the numerical convergence. In GRMHD models, we refine the AMR blocks as the central 
density grows to save the computational time. Typically the total number of AMR blocks increases 
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from ~ 0(H) 2 ) at initial to ~ O(10 3-4 ) at core bounce. We check whether this refinement procedure 
breaks such as the Hamiltonian constraint or the ADM mass conservation. 

In Fig J16l we plot the central lapse a, the error of the Hamiltonian constraint C% ( x 10) (see, Eq. 
(I54D ) and the ADM mass normalized by its initial v alue of magnetorotational collapse of a WOMq 
star with zero metallicity (lUmeda Nomotd 120081 ) . Such a highly massive star is considered to 
form BH and a good objection to test our GRMHD code. Initial conditions and the adopted EOS 
are the same as those used in model "GB12R020Sf" (see, SecJT]) and the minimum cell width is 
Ax = 600m. We also plot the deviation of the total angular momentum along the rotational axis 
from its initial value \AJ Z \ = \(J Z — J z ,o)/Jz,o\, here J z is defined by Eq. ([59]) and J Zt o is the value 
at t = ms. In this test, the total number of AMR blocks increases from 792 to 11432 through 
three times refinement procedures until the time of core bounce. The ADM mass and the total 
angular momentum should be almost constant until the core bounce and for a short while after it. 
This is because it takes 5000km/c ~ 10ms till the gravitational radiation, emitted at core bounce 
around the center, reaches the outer boundary 5000km. From this test, we find that the total 
angular momentum are conserved within ~ 1% until the time of core bounce. In addition, the time 
evolution of the central lapse is smooth and no influence of the refinement procedures is seen. As 
for the Hamiltonian constraint, C% is kept within several percentage until the time of core bounce 
except t = 28 ~ 31 ms and after a short while from the core bounce. We consider that the sudden 
increase of C% during t = 28 ~ 31 ms is due to the low resolution and thus the next refinement at 
t ~ 31 ms suppresses the error. After the core bounce at t ~ 32 ms, C% increases gradually due to 
the collapse bu t not rapidly and finally the calculation is crushed. We found the apparent horizon 
(|Shibatalll997l ) is formed at the end of the calculation and we thus consider the BH is born. 
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Fig. 16. — The Hamiltonian constraint C%(xl0) (green), the central lapse a (blue), the deviation 
of the total angular momentum |AJ 2 |(xl0) (red) and the ADM mass normalized by its initial value 
(black) are plotted against time. Refinement procedures and enforcing the Hamiltonian constraint 
are done four times at when arrows point in the top. Initial ADM mass is 4.337M and the total 
number of AMR blocks is increased as 792, 848, 3648, and 11432. 



In Fig. [T71 we display magnified views of Fig. [16] around the time of core bounce with two 
different numerical resolutions. We put the outer boundary at 5000 km for lower resolution model 
and 4000 km for higher one and maintain the AMR structure almost the same in both models. 
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Then the cell widths of the higher resolution model are 0.8 times smaller than those of lower 
resolution model. In Fig. \T7\ dashed and solid lines correspond to higher (minimum cell width is 
Ax = 480m) and lower (Ax = 600m) resolution model, respectively. From this figure, we see that 
the good numerical convergence is achieved. 
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Fig. 17. — Time evolutions of the central lapse (top-left), various total energies (top-right) and 
the maximum density (bottom-left) with different numerical resolutions. The minimum cell widths 
are set as Ax = 600m for solid lines and as Ax = 480m for dashed lines. The internal, kinetic and 
magnetic energies are normalized by 10 52 , 10 51 and 10 48 , respectively. 



7. Collapse of a 15M Q Star 

In this section, we calculate magneto rotational collapse of a 15M Q progenitor star as our 
practical test using GR/NMHD codes and check their abilities. We follow the gravitational collapse 
with varying the initial magnetic field, the stiffness of the EOS and the initial central angular 
velocity. We first describe our initial setups and then show results after subsection 17.41 



7.1. Equation of State 



We adopt a parametric type EOS (e.g.. iTakahara fc Satdll988l ) in this study. It is divided into 
two parts, "cold" and "thermal" part. The cold part is expressed as 



Pc 



K 2P r * 
K 3 p Ti 



< p < pi 
Pi< P< P2 

P2< P < P3 

P3 < P < Pa 



(76) 



where, we fix K\ = 2.78 x 10 14 in cgs units and -f^2,3,4 are determined from the continuity of P c at 
Pi, 2,3- P3 corresponds to the nuclear density p n uc- Polytropic indexes represent physical processes 
occurring during core collapse such as the electron-capture, onset of the neutrino-trap and the 
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nuclear repulsive force. We adopt two types of poly tropic indexes and, for convenience, we call 
"Sof t" and "Stif f ' EOS as summarized in Table [TJ Stiff EOS corresponds to the models reported 



m 



Mikami et al.l (|2008l ) and Soft EOS adopts smaller polytropic index compared to the Stiff case 



in the range of p\ < p < p2- That density region corresponds to the electron capture regime. 
The thermal part is expressed by 

P t = (r t - l)pe t 



(77) 



and we fix the index of the thermal part as T t = 1.3 in this study. Then the total pressure p and 
internal energy e contributed from both thermal and cold part are written by 



P = Pc + Pt 

e = e t +e c 

= £t + / ~17^ dp 
Jo P 



(78) 
(79) 



7.2. Grid Setup 

In NMHD models, we did not change their initial AMR structures and fixed them. On the other 
hand, we turn on the switch of AMR and refine AMR boxes in the vicinity of center in GRMHD 
models. This is because, the time step AT to keep the Courant-Friedrichs-Lewy condition (CFL 
condition) is determined from the maximum wave speed which usually becomes that of dynamical 
background and not hydrodynamical wave speed (i.e., the fast magneto sonic) in GRMHD models. 
The wave speed of the dynamical background is nearly the speed of light c and hardly change 
throughout the time evolution and also not depend on the hydrodynamical properties such as 
maximum density. We set AT as 0.3Ax m ; n /c in GRMHD models to keep the CFL condition where 
Ax m ; n is the minimum cell width in the computational domain. If we set maximum Lamr as 
the maximum allowed one (I/AMR,max = 8 or 9 in this study) from the beginning of calculation, 

Table 1. Parameters of cold part EOS 



i 




Soft EOS 






Stiff EOS 




Pi(g cm 3 ) 


Ki 


Ti 


Pi(g cm 3 ) 


Ki 


Ti 


1 


4 x 10 9 


2.78 x 10 14 


4/3 


4 x 10 9 


2.78 x 10 14 


4/3 


2 


1 x 10 12 


7.25 x 10 14 


1.29 


1 x 10 12 


4.66 x 10 14 


1.31 


3 


2 x 10 14 


3.17 x 10 14 


1.32 


2.8 x 10 14 


2.45 x 10 14 


4/3 


4 


1 x 10 16 


4.22 x 1(T 3 


2.5 


1 x 10 16 


3.42 x 1(T 3 


2.5 
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it takes too much time until corebounce since the time step depends solely on Ax m j n . Therefore 
in GRMHD models, we set maximum Lamr as 5 at the beginning and then increment it as the 
collapse proceeds to save computational time. We define the criterion to increment the maximum 
AMR level (i.e., refine the AMR boxes in the vicinity of center) that every time the central density 
exceeds 10 11 , 10 12 , 10 13 g cm -3 . Influences of box refinements are described in Sec. 16.4.31 

The outer boundary is taken at 5000 km from the origin and 8x8x8 AMR boxes with refine- 
ment level (L\mr = 0) cover the whole computational domain which is (x, y, z) = [—5000, 5000]km. 
We set the maximum allowed AMR refinement level to 8 in standard model and to 9 for high resolu- 
tion model, thus the highest resolution is Ax ~ 600m in standard and Ax ~ 300m in high resolution 
run. In standard models, the central region of (x, y, z) = [—60, 60]km is covered by Lamr = 7 and 
l 3 ^?/'- 2 ! ^5 30km is covered by Lamr = 8 (see, Fig lT8j) . On the other hand in high resolution run, 
the central regions of < 120, 60, 30 km are covered by Lamr = 7, 8, 9, respectively. Since 

the most dynamical and active region is above the surface of central core (r > 20km) and inside 
the prompt shock (r < 200km), we increase resolutions as much as possible in this region. Then 
the total numbers of AMR blocks become ~ 6000 in standard run and ~ 42000 in high resolution 
run. 




1 o' -5x10 6 o &xio s io 7 



Fig. 18. — Schematic figure of grid setup in standard resolution run. Each point represents center 
of cell. The central region of \x,y, z\ < 30 and 60 km is covered by Lamr = 8 and 7, respectively. 
In high resolution run, the central region of \x, y, z\ < 30, 60 and 120 km is covered by £amr = 9> 
8 and 7, respectively. 



7.3. Initial Setup 



Our progenitor is a 15Mq star with the metallicity Z=0.02 and the pre-collapse model is taken 
from lUmeda 8z Nomotol (|2008l ) which calculates stellar evolution with spherical symmetry. We add 
rotation and magnetic field to the spherical progenitor model. Since little is known about the 
rotational law and magnetic field configuration in the central iron core and its surroundings at the 
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pre-collapse stage, we assume the rotational law such as 



2 

n(w) = Q — p — (for NMHD) (80) 

Wq z + ZU 

u*ity = w 2 (n - n) (for GRMHD) (81) 



Where w = \J x 2 + y 2 , = y m 2 + u 2 and voq is the parameter which is set as 10 8 cm in this 
study. These configurations are commonly used rotational law in which the core, within ~ wq, 
rotates rigidly with angular velocity f^o and differentially beyond that. Both equations represent 
same rotational profile, since in the Newtonian limit of Eq. (|8ip becomes Eq. (|80p by replacing 
it* — > 1 and U0 — > 700^^ = tu 2 J7, where 7^ is a component of cylindrical flat metric. 

As for the initial magnetic field configuratio n, to ensure diverge nce-free constraint, we adopt 



the following form of vector potential (see, e.g., iTakiwaki et al.1 12009 ) . 



(A r ,A e ,A (j) ) = ( ,0,^-3^_ ro ) (82) 



Where Bq &z Rq are the parameters and r = \J x 2 + y 2 + z 2 . Rq is fixed as 10 8 cm in this study. This 
vector potential represents almost uniform magnetic field within ~ Rq and dipole-like magnetic field 
configuration beyond ~ Rq with the central magnetic field strength ~ Bq. We calculated several 
models with various initial magnetic field strength Bq, in Eq. ()82|) . and central angular velocity 
f2o, in Eq. (j80M8ip . We also add random perturbation to trigger asymmetric motion in the form of 
velocity with 5% amplitude when the maximum density exceeds 10 13 g cm -3 . Such perturbation is 
added in the central sphere of radius 10 8 cm via 

n = n ( 1+ 'ITw) (83) 

In this equation, Rq is a parameter fixed with Rq = 10 8 cm and —0.05 < a < 0.05 is a random 
number. Here, we have to comment about the momentum constraint, Eq. (|53p . in GRMHD model, 
since the momentum constraint would be violated by adding the perturbation. In our practical 
simulations, the momentum constraint Mi is not kept strictly during time evolution due to the 
numerical error. In our GRMHD model "GB12R020Sf", for instance, Cm x defined by Eq. ()84p is 
Cm x ~ 1-53 x 10~ 2 before the perturbation is added. 

C Mi = J P *Midx 3 (84) 

By adding the perturbation, this value increases to Cm x ~ 1.55 x 10 -2 . On the other hand, the error 
after core bounce is around Cm x ~ 4.5 x 10~ 2 in both models with and without the perturbation and 
we thus consider that the violation of the momentum constraint by the perturbation is negligibly 
small within the range of our numerical accuracy. 

Model names and adopted parameters are summarized in Table [2j The first character " N" 
and "G" of model names indicate that the calculation is done by NMHD and GRMHD code, 
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respectively. Numbers after "B" and "R" represent the exponent of magnetic field strength and 
the central angular velocity, respectively. The last characters St/Sf represents Stiff/Soft EOS 
adopted. As for the initial magnetic field strength, we adopted 0, 10 9 and 4.8 x 10 12 G so that we 
can easily check the roles of magnetic field. 

The initial magnetic field is first amplified mainly through the compression and the rotational 



wind i ng effects during colla pse from previous many studies (e.g.. lShibata et al.ll2006l ; lBurrows et al 



20071 ; iTakiwaki et al.l 1200911. For instanc e, the compression mechanism amplifies the magnetic field 
about ~ 10 3 times (jBurrows et alJl2007l l and the n B n ~ £>(10 12 1G is amplifi ed to ~ O(10 1 5 - 16 )G 
which is equivalent strength to that of magnetar ([Duncan &: Thompsonlll992l ). iHeger et al.l ((2005) 
studied stellar evolution with including magnetic field and rotation and they reported the strength 
of magnetic field is of the order of ~ 10(9)G or weaker and the poloidal magnetic field is much 
weaker, approximately 10 -4 , than the toroidal component. Thus our initial condition with purely 
poloidal and extremely strong magnetic field might be unrealistic one, however we employ such 
initial condition to see the effects of magnetic field easily and also to compare our results with 
other previous studies. In addition to the very strong initial magnetic field models, we calculated 
one model (NB09R02Sf) with initially weak magnetic field Bq = 10 9 G in a high resolution run. 
Since, there are possibly several non linear magnetic field amplification mechanisms in the vicinity 
of PNS such as MRI or dynamo mechanism which are intrinsically 3D phenomena, we examine how 
the initially weak magnetic field is amplified through this model. 



As for the central angular velocity fio we adopt 2 or 6 (rad/s). In lHirschi et al.1 ( 2004 1. they 
calculated the evolutions of various rotating stars with changing initial stellar mass and metallicity. 
Their results are that 25M@ star with solar metallicity (Z = 0.02) has Qq ~ 1.0 s _1 at the end 
of silicon burning stage. They also re ported that lowe r initia l metallicity raises the final angular 



veloci ty due to lower mass loss rate. lYoon Langerl (|2005l ) did similar works to IHirschi et al 



(12004! ) . though they included magnetic effects. Their results showed magnetic torques lower the 
local specific angular momentum approximately one order of magnitude compared to non-magnetic 
field cases. However these two works are done by one dimensional calculations and are not still 
conclusive results. Thus, even though our adopted parameters are comparatively faster than the 
theoretical works those values may be reasonable. 



7.3.1. Initial Setup for GRMHD 

In GRMHD calculations, we have additional setups to be done which are constraining the 
Hamiltonian and momentum equations (|52p and f)53[) . To solve Eqs. (|52M53p . we assume conformally 
flat metric for initial c ondition, in which 7^ = 5ij (therefore, F, = 0) and K = 0. Then, following 
(jShibata &: Urvul 120021 1 . constraint equations become the following 4 equations to obtain rest of 
BSSN variables cj> and Ay and gauge variables a and 



A 



-2-irSoe 



5«< 



(85) 
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A flat (ae*) = 27rae^(5 e- 6 ^ + 2 7 ^) + ^——6*6*1^ (86) 

o 

<%A flat /^ + ^\ ki - 2A ik 6 k i (aj - ^eO = WTraS^ (87) 

h = ^ U k /3 k d + 6 jk p k 7i - (88) 

Where Afj at is the Laplacian in flat space. We solve above 4 equations by iterative method with 
initial guess as <f> = 0, Aij = 0, a = 1, j3 % = 0. Then first, we evaluate conservative variables Q 
from given metrics and primitive variables P from a progenitor model. Second, we solve above 
each equation. We iterate these two procedures until sufficient convergence is achieved. After these 
procedures, the Hamiltonian constraint C% at initial is kept below arbitrary chosen small number 
(in our calculations 10~ 8 is adopted). 



7.4. Results 

7. 4-1- Global Dynamics 

We first show the time evolutions of the maximum density /0 max m Fig ]19l which tells us rough 
overview of how the core collapse proceeds. In this study our calculations are done in full three- 
dimension and the central density does not always become the maximum one, thus we do not use 
central value. However all of our models show that the maximum density exists nearly the center. 
In Fig ll9l each line corresponds to different model and the model names are shown in the bottom 
part. From this figure, we can find /0 max is increased ~ 30% in GRMHD models when compared 
with corresponding NMHD models (e.g., "GB12R020SF vs "NB12R020SP). We however do not see 
any rapid increase of p max which possibly indicates collapse toward BH formation. Another feature 
is that the central density gradually increases after core bounce when the magnetic field exists, 
on the other hand the central density decreases in non-magnetized models (e.g., black vs red solid 
curves). This difference means that the magnetic torque works strongly and the maximum density 




time („) 



NB00R020Sf NB12R060Sf NB00R020St GB00R020Sf 

NB12R020Sf NB09R020Sf ,NB12R020St . GB12R020Sf 

Fig. 19. — Time evolutions of the maximum density in different models. Model names are listed in 
the bottom. 

increases due to the angular momentum transfer. In rapidly rotating model "NB12R060Sf ' , the 
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maximum density at core bounce marginally exceeds the nuclear density (p nuc = 2 x 10 14 g cm -3 ) 
and it is thus not rotational supported core bounce but due to the nuclear repulsive force which 
is the same as the rest of models. In this model, p max eventually relaxes to similar value to other 
models after several oscillations. We also summarize physical properties at core bounce in Table 

El 

Fig. [20] displays time evolution of density profiles along the x axis in model NB12R020Sf(/e/t) 
and NB12R020St(ri(7/i£). Numbers beside each line denote time in ms. Prompt shock which is 
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Fig. 20.— Density profiles along the x axis in model NB12R020Sf(Ze/i) and NB12R020St(r^/it). 
Numbers denote time in ms. 

formed at core bounce moves outward and then stays around x ~ 200-300km in model NB12R020Sf, 
on the other hand the shock moves further out in stiff EOS model NB12R020St. Model NB12R020St 



adopts same EOS as that of iMikami et alJ (12008I ) and our result seen in the prompt shock propa- 



gation agrees well with their results. 

In Fig J21l we display time evolutions of the rotational, internal and magnetic energies in left 
four panels and we also plot comparison between GRMHD(GB12R020Sf ) and NMHD (NB 12R020Sf ) 
models in right two panels. The magnetic energy of model "NB09R020SP is too small compared 
to other strong field models and we display it separately with different range in bottom-right panel. 
We see that the rotational and the internal energies are kept almost constant or gradual increase 
after core bounce, on the other hand the magnetic energy increases rapidly, approximately ~ 2 
orders, after core bounce. In strongly magnetized models (upper-middle panel), the final magnetic 
energies saturate around ~ 5 x 10 50 ergs. When we compare GRMHD and NMHD models, shown 
in right two panels, the evolution tracks look similar except E rot and E mag after t > 80ms. We 
consider the difference is originated from the bipolar-outflow and will be described in Sec. I7.4.2i 



Finally we compare our results shown here with those reported by other groups. lObergaulinger et al 



(|2006l ) reported magnetorotational collapse in axisymmetry with various initial rotation, magnetic 
field and EOS and also with including general relativistic effects by replacing spherical Newtonian 
potential with "Tolman-Oppenheimer-Volkoff" potential. They showed that maximum rest mass 
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Table 2. Initial parameters and adopted EOS 



model name 


So(G) 


Ho(rad/s) 


Anag 




Egrv (ergs) 


a c 


M bar /M 


Madm/Mq 


EOS 




NMHD 


NB00R02St 


0E0 


2 


0E0 


0.135 


-6.46E51 




2.0483 




Stiff 


NB12R02St 


4.8E12 


2 


2.95E-4 


0.135 


-6.46E51 




2.0483 




Stiff 


NB00R02Sf 


0E0 


2 


0E0 


0.135 


-6.46E51 




2.0483 




Soft 


NB12R02Sf 


4.8E12 


2 


2.95E-4 


0.135 


-6.46E51 




2.0483 




Soft 


NB12R06Sf 


4.8E12 


6 


2.95E-4 


1.22 


-6.46E51 




2.0483 




Soft 


NB09R02Sf 1 


1E9 


2 


1.15E-11 


0.135 


-6.46E51 




2.0483 




Soft 




GRMHD 


GB00R02Sf 


0E0 


2 


0E0 


0.132 


-6.62E51 


0.994 


2.0688 


2.0673 


Soft 


GB12R02Sf 


4.8E12 


2 


2.00E-4 


0.132 


-6.62E51 


0.994 


2.0688 


2.0673 


Soft 



This model is calculated with high resolution. 



Note. — Each column denotes model name and initial parameters. From left; model name; magnetic field Bo; central 
angular velocity Qo; /3mag = -E mag / 1 E gTV | ; /3 r ot = T ro t / 1 E gTV | ; gravitational energy E grv ; central lapse; baryon rest mass 
M^/Mq; ADM mass M ADM /M ; and adopted EOS. 



Table 3. Physical properties at core bounce 



model name 


Pmax/10 14 (g cm- 3 ) 


/^m ag 


Aot(%) 






NMHD 


NB00R02St 


3.54 





3.62 




NB12R02St 


3.55 


6.17E-4 


3.58 




NB00R02Sf 


3.95 





3.00 




NB12R02Sf 


3.96 


5.88E-4 


2.82 




NB12R06Sf 


2.25 


165E-1 


13.9 




NB09R02Sf 


3.98 


2.88E-11 


3.04 




GRMHD 


GB00R02Sf 


4.87 





2.25 


0.839 


GB12R02Sf 


4.87 


2.31E-4 


2.20 


0.839 



Note. — From left; maximum density p m ax; Anag = -Emag/|-Egrv|; 
A-ot = T r ot / 1 E grv | ; central lapse a c at core bounce. 
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Fig. 21. — Left four panels; Evolutions of rotational, internal and magnetic energy. Initially 
weak magnetic field model "NB09R020Sf is plotted with different range in the lower middle panel. 
Right two panels; Comparison between GRMHD(GB12R020Sf) and NMHD(NB12R020Sf) models. 
Upper and lower panels represent the maximum density and each energy component, respectively. 



density is increased several 10 % after core bounce when they compare GR and Newtonian models 
and also that the magnetic field works to raise the maximum rest mass density. These features agree 
to ours since ~ 30% rise in the maximum density in our GRMHD model can be seen. Additionally, 



i f we compare our results with previous three-dimensional NMHD work reported by iMikami et al 



(|2008l ). similar time evolutions are obtained such as time evolution of various energy components 
and also the shock propagation (as seen in Fig. [20|h Thus, we consider the results shown here are 
common and robust features. 



7.4-2. Formation of Outflow 

Next, we describe the formation of bipolar outflow. In all of our strongly magnetized mod- 
els, bipolar outflow is formed in a similar manner and we thus present mainly one representative 
model "NB12R020Sf" in this subsection. Figj22l and Figj23l show the density contour in model 
"NB12R020SP at different time slices. FigJMl and Fig US are the same as Figj22] but are with 
the color coded contour of plasma beta (/3 P = P ga s/Pmag) in logarithmic scale and the flow ve- 
locity in white arrows. Black curves represent the iso-density contour. The depicted region is 
(x,y,z) = [— 150, 150]km and, in each panel, bottom-left, top-left and top- right part represents 
xy (equatorial), xz and yz plane, respectively. As shown in these figures, the strongly magne- 
tized regions where log /3 p reaches ~ appear along the rotational axis, which means the magnetic 
pressure is comparable to the matter pressure. Then high velocity outflow is launched along the 
rotational axis (see, Figl26l for NMHD and Figl27l for GRMHD) while inflow appears along the 
equatorial plane. On the other hand, non-magnetized models and weakly magnetized model 
"NB09R020Sf" do not form any bipolar outflow as seen from two right panels of Figl261 
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Fig. 22. — Logarithmic scale of rest mass density in model NB12R020Sf at different time slices are 
depicted. Time slices are chosen at t=66ms (nearly the time of core bounce) and t=74ms. 
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Fig. 23.— Same as Fig. [22] but for t=82ms and t=90ms. 
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Fig. 24. — Logarithmic scale of plasma beta in model NB12R020Sf at different time slices are 
depicted. Time slices are chosen at t=66ms (nearly the time of core bounce) and t=74ms. Arrows 
represent flow velocity and black curves represent iso-density contour. In the figure, we cut-off 
log/3 p and the flow velocity higher than 3 and 3 x 10 9 cm s _1 , respectively. 
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Fig. 25. — Same as Fig. [Ml but for t=82ms and t=86ms. 
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Fig. 26. — Three dimensional pictures of bipolar-outflow. Left two panels are "NB12R020SP 
and right two are "NB00R020Sf" and \x,y, z\ < 300km is drawn. The central spherical-like shells 
represent iso-density surfaces and the outermost shell (green) corresponds to 10 9 g cm -3 . White 
lines are the magnetic field lines. Translucent and opaque red surfaces are iso-radialvelocity surfaces 
and each corresponds to v T = 10 9 cm s _1 (translucent) and v r = 2 x 10 9 cm s _1 (opaque). In non- 
magnetized model "NB00R020Sf" , high velocity bipolar outflow did not appear. 





Fig. 27.— Same as Figj26] but of GRMHD models (left; "GB12R020SP and right; 
"GB00R020Sf"). Both snapshots are taken at t=87 ms. 
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When we compare GRMHD(GB12R020Sf) and NMHD(NB12R020Sf) models, we see that the 
shock front of bipolar outflow moves faster, approximately a factor of 2, in GRMHD model as 
can be seen in Fig. [28l Since the bipolar outflow is driven by the angular momentum transfer 
(described in the next section), higher velocity outflow reflects that GB12R020Sf extracts larger 
angular momentum compared to NB12R020Sf. Then the magnetic energy is increased while the 
rotational energy is decreased as seen in right-bottom panel of Fig. [2TJ 



Fig. 28. — Time evolutions of the bipolar shock front along the rotational axis (Only those of 
north hemisphere are shown). Solid and dashed lines are results of NB12R020Sf and GB12R020Sf, 
respectively. 



7.4-3. Driving Mechanisms of Outflow 

In this section, we describe the driving mechanisms of the outflow. As mentioned in previous 
subsection, all strongly magnetized models exhibit high velocity outflow along the rotational axis. 
The ultimate energy source of outflow is the angular momentum transfer from the central object 
which can be seen from Figl291 In this figure, the angular velocities along the x axis of different 
models at different time slices are shown. For instance, from left two panels, we can see the angular 
velocity within x < 10km decreases rapidly, a factor of several, from t = 66ms to t ~ 80-90ms. 
On the other hand, non-magnetized models (right two panels) do not show any deceleration (note 
that kink in the angular velocity profile at x ~ 3km shown in t = 120ms of "NB00R020Sf is 
originated from such as meridional circulation or displacement of the mass center). From this fact, 
the angular momentum is extracted from the central object by the magnetic field. The extracted 
angular momentum is first converted to the magnetic field mainly via the magnetic field lapping. 
Then, there are two types of driving mechanisms of the outflow, the magneto-spring and the 
magneto-centrifugally supported mechanisms. We consider, from left two panels of Fig l26l that 
initial mechanism is the magneto-spring effect and then it transitions to the magneto-centrifugally 
supported mechanism. This is because from FigJ26|, we see the magnetic field lines are highly 
twisted inside the shell at the onset of launching the outflow (t=120ms), however this twisted 
configuration is stretched eventually, as if the compressed spring would do, toward the north-south 
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Fig. 29. — The angular velocity profiles along x axis at different time. From top-left to clock wise 
direction, "NB12R020Sf", "NB00R020SP , "GB00R020SP and "GB12R020Sf" are displayed. Note 
that in all four models shown in this figure, the time of core bounce is approximately 66ms. Kinks, 
e.g., shown around x ~ 3km in the top-left panel at T = 120ms, represent the retrograde of the 
angular velocity. 



direction as seen in t=171ms panel. Final magnetic field configuration is less toroidal dominant 
compared to that of t=120ms and the matters stream away along these helical magnetic field lines. 
This is also seen in Figl30l which shows time evolutions of the outflow velocity (V z ), the toroidal 
magnetic field (B^) and angle 0=cos -1 ((B • ~V)/BV) between the magnetic field and the velocity 
vector, along the rotational axis of model "NB12R020Sf". From upper-right panel, we find the 
strongly amplified B^ at t ~ 120ms and z < 80km which shortly disappears. At the same time, the 
high velocity region appears in the upper-left panel. The angle 9 is nearly ~ 90° at z < 50km and 
t > 120ms which indicates that the outflow is driven by the gradient of magnetic pressure (i.e., the 
magneto spring effect). Then V z is accelerated up to z ~ 110km along the magnetic field (9 ~ 0°, 
i.e., the magneto centrifugal effect). 



In 



Shibata et al.l (|2006l ). they calculated 2D axisymmetric GRMHD simulations and reported 
MHD outflow is first driven by the magneto-spring effect and eventually by the magneto-centrifugally 
supported mechanism. Thus our results of launching processes of MHD outflow are qualitatively 
similar to theirs. 
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Fig. 30. — Time evolutions of the velocity component (V z , upper-left), the toroidal magnetic field 
(Bff,, upper-right) and angle 6 (bottom) between the magnetic field and the velocity vector, along 
the rotational axis of model "NB12R020Sf ' . 



7.5. Non-Axisymmetric Motion 

As described in the previous subsection, our results seen in dynamical evolutions are qual- 
itati vely the same as tho s e reported in previ o us 2D axisymmetric M HD works (for NMHD see, 



Kotake et all (120051): ISawai et all (120051): burrows et al.1 (|2007l ) and for GRMHD see, e.g., 



Obergaulinger et al.l (|2006l ): IShibata et al.l (|2006l )) in the way that the equatorial inflow and the 
bipolar outflow along the rotational axis due to the magnetic field. We also calculated several mod- 
els with tilted magnetic field axis against the rotational axis to induce larger non-ax isymmetry and 



found the outflow is driven along the rotational axis similar to lMikami et al.l (120081) . According to 



many previous studies, both in GR and the Newtonian limit (e.g.. lOtt et alJl20071 ; IScheidegger et al 



2009), the nascent neutronstar is sensitive to the rotational instability predominated by the "m=l" 
non-axisymmetric mode within our initial rotational parameter range. To confirm our code can 
actually reproduce some nonaxisymmetric modes characteristic to rotational collapse of m assive 
stars, we monitor the non-axisymmetry with the same approach as IScheidegger et al.1 (|2009i ). We 
decompose the density into the Fourier components along the equatorial ring of radius 40km which 
is is beyond the central rigidly rotating region and inside the prompt shock (see, Fig|20p. Fourier 
amplitude of mode "m" is defined as the following equations. 



p(m,z,(f>) 



E A m ( 



zd, z)e 



im<f> 



m= — oo 



1 



A m (w,z) = — I p(zD,z,(j))e 



2tt 



(89) 
(90) 
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In Figl31[ we plot the time evolutions of normalized mode amplitude |A m |/|Ao|. Prom this figure, we 
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Fig. 31. — Time evolutions of the normalized mode amplitude |A m |/|Ao|. Model names are listed 
in the bottom part. 



see that m=4 mode is the most dominant mode before core bounce (t ~ 66ms) since our cartesian 
grid induces quadrupole numerical noise at initial. However, the linear amplification phase starts 
immediately after core bounce and several ms later it reaches the non-linear phase (t > 75ms). 
During the non- linear phase, dominant mode becomes m=l and their normalized amplitu de exceed 



x", 0.1. This is consistent with the structure, a so-called one-armed spiral structure seen in lOtt et al 
(|2005l ). Since j3 m t at core bounce reaches ~ 3% and also keeps > 1% after core bounce, we consider 
that the low-|T/W| instability causes this non-axisymmetric configuration. Like these, even though 
the outflow structure is almost axisymmetric from broader perspective, non-axisymmetry develops 
and show significantly large mode amplitude in the vicinity of center in the self gravita ting system. 
This non-axisymmetry may alter the gravitational wave form (jScheidegger et al.l l2009). 

Next, we describe about the amplification of initially weak magnetic field (Bq = 10 9 G) in 
model "NB09R020SP. As mentioned above in Sec l7.4.2"l only through the field- wrapping and 
the compression mechanism, the magnetic field cannot be amplified strongly enough to drive the 
outflow soon after the core bounce as seen in other models with initially strong magnetic field. 
However, there might be several magnetic field amp lification mechanisms to be occurred after core 



bounce such as the MRI o r the dynamo mechanism (jAkivama et al. 



2003 



Cerda-Duran et al. 



2007 



Obergaulinger et al.ll2009l ) in addition to the aforementioned linear mechanisms. If some of these 



mechanisms operate within dynamical time scale, the saturated magnetic field is considered to 
possess enough capability to affect the explosion dynamics. The largest difference between 3D and 
2D (axisymmetric) in the amplification of the magnetic field is conversion from toroidal to poloidal 
magnetic field, since the toroidal to poloidal conversion can never happen in 2D axisymmetric 
motion. This is because, from the Faraday's law, time evolution of a poloidal component B po1 
becomes 



d t B po1 + d pol± (B pol v pol± -v pol B pol± ) 
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+ d tor (B pol v tOT -v pol B toI ) = (91) 

Here, "pol" and "tor" represent a poloidal and toroidal component, respectively, and "pol_L" is a 
perpendicular one to both the "pol" and "tor" components. In axisymmetry, 9 t or = and thus 
B tor cannot be converted to B po1 , however in full 3D, several non-axisymmetric fluid motions (e.g., 
the Parker or the Tayler or the convective instabilities) let dt or 7^ and close the conversion cycle 
(i.e., from poloidal to toroidal and toroidal to poloidal). Through our weakly magnetized model 
"NB09R020Sf", we examine how the magnetic field is amplified after core bounce and whether the 
amplified magnetic field affect the explosion dynamics or not. 

In Fig J32l we display time evolution of the magnetic field strength in logarithmic scale of 
model "NB09R020Sf" which is high resolution run. Within r < 20km, strength of the magnetic 
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Fig. 32. — Time evolution of the magnetic field in logarithmic scale of model "NB09R020Sf" . 

field shows stratified configuration compared to r > 20km and is amplified strongest (~ 2 x 10 13 G) 
among the numerical domain. Since, within 10 < r < 20km, matters rotate differentially with the 
steepest angular velocity gradient (see, Fig|29p. the magnetic field strength is higher than the value 
of central region r < 10km where the rotation is almost rigid. Magnitude of the amplification is 
~ 10 13 /10 9 = 10 which is close to the value predicted by compression and we thus consider the 
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dominant field amplification mechanisms within 10 < r < 20km are compression and the rotational 
winding effect. 

In contrast, we see the region r > 20km is highly non- homologous and the magnetic field 
configuration pattern changes momentarily, however the maximum strength keeps the same level 
of the order of < 10 13 G throughout our calculation (~ 40ms after core bounce). The stochastic 
configuration pattern beyond r ~ 20km is mainly triggered by the entropy driven convection. The 
reason is like this. There are several candidates to configure such flow pattern such as the convective 
motion, the MRI, the Parker instability (or the magnetic buoyancy) or the Tayler instability. 
However, among them, the growth time scales of the Tayler and the Parker instabilities are too 
long and they are inefficient during our calculation time. For instance, the growth time scale of the 
Tayler instability is the or der of the Alfven crossing time of the system and ~ 0(1) h according to 
Cerda-Duran et al. (|2007l ) in a weakly magnetized limit. As for the Parker instability, the growth 



time scale (r mag ) can be esti mated by using the frequency of the magnetic buoyancy N ma , g as 



2nN m } lg . According to lAchesonl (|1979l ). iV mag is defined by using the the Alfven velocity 



c„4, the speed of sound a, the toroidal magnetic field B& and the rest mass density p as 



A^a g = ^g-v(hA) (92) 



Where, "y=dlnp/dliip\s, S is the entropy and g is the gravitational acceleration. In what follows, 
we adopt a pseudo-entropy defined by S=P t /P c - In Fig. [33J, we display color coded contour of 
(right) in addition to the Brunt- Vaisala frequency TV 2 (left) defined by 



dS 



■ VS (93) 



in log scale of model "NB09R020Sf ' . In each panel, color- less area is where N 2 /N^ has positive 
value and is thus stable region against each mechanism. The specific angular momentum with 
positive gradient with respect to r has stabilizing effect on the convective motion, however it 
has negative gradient in our models and the unstable region becomes larger when we consider 
the contribution from it (i.e., the Solberg-H0iland criterion). From Fig. [33j we see the PNS is 
convectively and magnetic buoyantly unstable. However, the growth time scale of each mechanism 
differs widely, r ~ 27riV _1 ~ 1 — 10ms for convection and T mag ~ 27rAT-i g > 0.1 - 10s for magnetic 
buoyancy. Since our calculation time is ~ 40ms after core bounce, the Parker instability does not 
grow while the convection can grow sufficiently within our simulation times. The convective-dynamo 
is thus considered to contribute to the magnetic field amplification mechanism. 

We also examine the possibility of the MRI. Since, from local linear analys is, the MRI would 



occur when the rotation is differential with negative angular velocity gradient (jBalbus &: Hawlev 



199ll ). From Figj29j we see inside the shock (r < 100km) has negative angular velocity gradi- 
ent and is thus unstable against the MRI. However, to follow the MRI by numerical simulation, 
the critical wave length of the MRI "Amri" has to be resolved at least ~ 10 numerical grids 
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Fig. 33. — Logarithmic scale of the Brunt- Vaisala frequency — iV 2 {left) and the magnetic buoyant 
frequency -iV 2 ag (right) of model "NB09R020SP. Color-less area is where N 2 /iV^ ag has positive 
value (i.e., stable region) and white curves are the iso-density contour. 



(jObergaulinger et al.1 120091 ) . Here, Amri is defined by 



Amri = 47tca -m 



on 2 

dm 



2\ -1/2 



(94) 
(95) 



In FigJ34|, we display the ratio of the critical wave length Amri to the local numerical grid width 
in left panel and rough estimation of the growth time scale of the fastest growing MRI mode 
"TMRi(ms)" in log scale in right panel, tmri is defined by (see. iBalbus fc Hawlevlll99ll ). 



r MRI 
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m- 



dn 
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dm 



-i 



(96) 



From left panel of this figure, we see the most part has negative value and therefore the MRI 
cannot be resolved. We have to employ > 10 times higher resolution (Ax < 30 — 60m) to resolve 
the MRI or, if we adopt ~ 10 times larger initial magnetic field (~ 10 10 G), we possibly manage 
to do it with our high resolution run marginally, since the wave length Amri is proportional to 
\B\. Here, we comment about our strongly magnetized models (Bq ~ 10 12 G). If we extend what 
we mentioned just above (i.e., to adopt stronger initial magnetic field), we may easily resolve the 
MRI since the critical wave length Amri is approximately 10 3 times larger compared to weakly 
magnetized model. From Fig. [34"1 we can estimate the wave length as Amri ~ 60m (yellow region) 
and Amri ~ 600m (orange region) for 30 < r < 60km where the resolution is Ax ~ 600m. Then 
in the strongly magnetized models, Amri multiplied 10 3 becomes Amri ~ 60 — 600km. However, 
the system scale (i.e., inside the prompt shock) is ~ 200km and therefore those modes which have 
larger wave length than the system cannot last. In case an MRI mode lasts and if we can resolve 
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Fig. 34.— Contour of Amri/Ax (Ze/f) and t M ri (ms) (ri^fa) in log scale from model "NB09R020Sf ' . 
Colorless areas in left panel are where the MRI stable regions (i.e., positive angular velocity 
gradient). 



it, the magnetic field soon reaches saturation strength in several times of the growth time scale. 
Akivama et al.l (J2003J) derived the saturated magnetic field strength £? sa t,MRi as 



-^sat.MRI ~ 



(97) 



and it becomes ~ 10 15 G with our initial rotational parameter. All of our strongly magnetized mod- 
els exhibit saturated magnetic field strength of the order of ~ 10 15 ~ 16 G after c ore bounce (the valu e 
is consistent with t hose reported by prev ious many stu dies, e.g., for NMHD iKotake et al.l (120 05): 



Sawai et all (12005TI : Burrows et all (12007 ). for GRMHD. lObergaulineer et all ([20061 ) : 



Shibata et al 



(|2006l ) and for 3D works. iMikami et al.l (|2008l ) ; IScheidegger et al.l (|2009l )). The value is comparable 
to -B sati MRi and we thus consider the initially strong magnetic field is first amplified by the com- 
pression and the winding effect with the amplification magnitude of the order of ~ 10 3 and then the 



MRI operates to amplify the magnetic field up to the saturation strength (in lObergaulinger et al 



2006 



Shibata et al.ll200q . they reported the MRI operates with adopting similar initial magnetic 



field strength ~ 10 12 G). However, since just only the linear amplification mechanisms amplify the 
magnetic field up to ~ 10 15 G which is close to the MRI saturation level, to see the effects of the MRI 
amplification more clearly, we have to adopt sufficiently weak magnetic field (e.g., -Bo ~ 10 10 G) 
which does not reach -B sa t,MRi only through the linear amplification mechanisms but sufficiently 
strong to resolve by numerical simulation. 

As for the calculation time, the length of 40 ms after core bounce is marginally sufficient for 
the inner region (r < 40km) from the right panel of Fig. 134"! however beyond that region we have 
to evolve more than several hundreds ms. If we capture the linear amplification, it may reaches 
the saturation phase within the several rotational periods. At this saturation phase, whether the 
magnetic field is sufficiently strong to affect the explosion dynamics and how the magnetic field 
configuration is cannot be clarified without numerical simulation and this would be our future work. 
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8. Summary and Discussions 

The explosion mechanisms of the core-collapse supernovae have been unknown and fascinating 
problems for several decades. Recent observations show several common features seen in the CCSNe 
that some types of them are bipolar like and sometimes non-axisymmetric explosions. Motivated 
by these, we now have to take into account the effects of asymmetry into numerical works to 
uncover the explosion dynamics. Fortunately, recent development of computational resources enable 
us to handle the numerical simulations in the context of three dimension. We therefore have 
developed two types of three dimensional magneto hydrodynamical codes. One is in the Newtonian 
approximation (NMHD) and the other is in the full general relativity (GRMHD). The features 
of our codes are; (1) Adoptive Mesh Refinement to cover the wide dynamical ranges; (2) high 
resolution shock capturing scheme with Roe-like (in NMHD) and HLL (in GRMHD) flux; (3) 
several reconstruction schemes to maintain high spatial resolution; (4) time update of the matters 
and the metric is done by the iterative Crank-Nicholson scheme; (5) the constrained transport 
to evolve the magnetic field; (6) any types of the EOS can be adopted; (7) the Poisson solver 
with BiConjugate Gradient Stabilized Method under our AMR structure to solve the self gravity 
(in NMHD) and the non-linear Poisson like equations for the Hamiltonian and the momentum 
constrains (in GRMHD). 

In this paper, we described our numerical methods in detail and did several tests to con- 
firm their abilities through the simple shock tube tests; the Poisson solver for the spherically 
distributed matters; conservation of the mass, the energy and the local/global angular momentum; 
the quadrupole linearized gravitational wave; the rotating neutronstar in equilibrium states and so 
on. Through these tests, we confirmed that our codes reproduce the numerical error convergence 
predicted by our adopted reconstruction schemes and also confirmed that the accuracy of our code 
is sufficient to follow the dynamical evolution of CCSNe. And as for the first test of CCSN simula- 
tion, we calculated collapse of a 15M Q progenitor with varying the initial magnetic field strength, 
the angular velocity and the stiffness of the Polytropic EOS by our GRMHD and NMHD codes. 
Our main results and some discussions are as the following. 

(1) After a short while (~ 20ms) from the time of core bounce, high velocity (V r ~2x 10 9 cm 
s _1 ) bipolar outflow is driven from surface of the proto-neutronstar (\z\ ~ 30km) and moves through 
along the rotational axis. The bipolar outflow does not appear in the non-magnetized and the 
initially weak magnetized models which indicate the outflow is magnetically driven outflow. The 
energy source of this outflow is the extracted angular momentum of the central proto-neutronstar 
which is transfered by the magnetic torque. The driving mechanisms are first by the magneto 
spring effect and then we consider the magneto centrifugally supported outflow. 

(2) In our self gravitating system, the non-axisymmetry develops immediately after the core 
bounce with the linear amplification at first and soon reaches the non-linear phase. The dominant 
non-axisymmetric mode is the m=l mode during the non-linear phase and the one-armed spiral 
structure is also can be seen. Since our initial rotational velocity (Jl c > 2 rad/s) satisfies the low- 
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|T/W| instability criterion (/3 > 1%) after core bounce, we consider that these non-axisymmetric, 
spiral mode is originated from the low-|T/W| instability. However, these non-axisymmetric motions 
are confined in the vicinity of the center and, in general terms, the global structure of bipolar outflow 
is qualitatively the same as those reported in previous 2D axisymmetric MHD works in the way 
that the equatorial inflow and the bipolar outflow along the rotational axis. 

(3) In weakly magnetized model in which the initial central, poloidal magnetic field is 10 9 G, the 
convective over turn highly deforms the magnetic field configuration. However, with our resolution 
of Ax ~ 300 — 600m and limited computational time (~ 40ms after core bounce), we did not 
find the exponential growth of the magnetic field which can be seen if the magneto-rotational 
instability works. If we employ 10 times higher resolution or 10 times stronger initial magnetic 
field, we possibly capture the MRI marginally. Time scale of the MRI is tmri ~ 0(1) ms inside 
r < 30km which is comparable to the dynamical time scale and is sufficiently short to follow by the 
numerical simulation. Even if the MRI operates in a weakly magnetized model, whether the MRI 
amplifies the magnetic field strongly enough to affect the explosion dynamics or not and whether 
the amplified magnetic field contributes to launch the outflow are big issues. Especially, since the 
the amplified magnetic field through the MRI would be less directional (i.e., the magnetic field is 
not amplified intensively along the rotational axis as seen in Fig|26l & |27|) . we have to examine 
whether the amplified magnetic field affect the mangeto-rotational explosion scenario. 

(4) By comparing GRMHD and NMHD models, we found that the gravitational effect works 
a little bit stronger in GRMHD models which can be seen in the increase of the central density as 
~ 30%. However, the global dynamical evolutions are similar such as the time of core bounce and 
formation of the bipolar outflow. Therefore we consider that the Newtonian approximation for the 
low mass range (< 15M ) is acceptable at least for several ten ms after core bounce. In this study, 
our progenitor is a 15M Q star which is small among the mass range of CCSNe progenitors and 
thus, if we adopt much larger mass such as ~ 40 — lOOM©, the general relativistic effects become 
stronger and the qualitative differences may be appeared even in a similar time scale as those used 
in this report. Such high mass range calculations are now in progress and will be reported in near 
future. 

To confirm the validation of numerical results, we have to connect them to the observations. 
One important observed object is the gravitational wave. Since we cannot observe directly in 
the vicinity of the center of CCSNe by the electro-magetic wave, the gravitational wave is one of 
the limited ways with which we can observe directly. In this report, though we do not evaluate 
the gravitational wave forms, the non-axisymmetric motion and the bipolar configura ti on ap pear. 



These motions alter the gravitational wave forms as reported by IScheidegger et al.l (|2009l ) and 
we will examine the effects of, e.g., the progenitor mass or the non-axisymmetric motions or the 
magnetic field in the context of full general relativity. This will be our future work. Another object 
to confirm is the ejected elements accompanied with the explosions. CCSNe eject abundant heavy 
elements which are synthesized during the progenitor's main s equence age to the i r final fate and 



also during the explosion via, such as, r-process nucleosynthesis (jWanaio et al. 2002; F uiimoto et al 
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20071 ; iKuroda et alj|2008l ). However, the ejected chemical compositions and abundances depend on 



the detailed properties of the explosion dynamics and we can still not explain the observed chemical 
abundances, one reason is due to the lack of comprehension about the explosion dynamics. By 
comparing our numerical results with the observations, we can feed back the observational studies 
to our numerical models and input physics. 



Numerical computations were carried out on Cray XT4 at Center for Computational Astro- 
physics, CfCA, of National Astronomical Observatory of Japan. This work was partly supported 
by Grants-in-Aid for JSPS Fellows and for the Scientific Research from the Ministry of Education, 
Science and Culture of Japan (20105004). We are grateful to an anonymous referee for his/her 
valuable and constructive comments. 
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